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ABSTBACT 


This  thesis  treats  the  problem  of  determining  the  eleva- 
tion angle  needed  to  initialize  an  underwater  sound  ray 
tracing  algorithm  used  to  locate  the  position  jf  a  target 
vehicle.  At  regularly  spaced  time  intervals  the  vehicle 
pings  a  synchronized  sound  signal  wnich  is  received  by  a 
(short  base  line)  sonar  array  containing  four  hydrophones 
positioned  at  four  of  the  corners  of  a  cube.  Tie  wavefrcnt 
direction  angles  are  determined  from  the  arrival  times  at 
the    four   hydrophones. 

Current  methods  for  using  such  time  data  t:>  produce  an 
apparent  position  suitable  for  ray  tracing  are  reviewed. 
Then  four  new  methods  are  developed  and  documented  mathemat- 
ically. All  methods  are  compared  under  a  simulated  environ- 
ment of  a  sound  speed  profile  which  is  linear  with  depth. 
One  of  the  new  methods  is  judged  to  be  an  impr>  vement  over 
current  methods  in  this  idealized  environment.  Finally  the 
improved  method  is  used  to  estimate  the  variability  in  the 
time    data   from    a  real  hydrophonic    tracking    problem. 
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I«    IHTBQIMJCTIQN    AND    BACKGROUND 

A.  BACKGROUND 

Suppose  that  an  underwater  acoustical  source  emits  a 
signal  at  a  known  time.  Suppose  that  the  times  of  arrival 
of  that  signal  at  the  hydrophones  on  a  three  dimensional 
sensing  array  are  also  known.  Finally  suppose  that  the 
speed  of  scund  in  water  is  modelled  as  being  homogeneous 
over  time  and  horizontal  displacement,  varying  only  with 
depth.  Then  the  angle  of  elevation  (A),  and  the  time  (T)  of 
the  arrival  of  the  signal  at  the  acoustic  center  of  the 
hydrophone  array  can  be  determined.  If  the  relationship 
between  the  speed  of  sound  and  the  depth  under  water  is 
known  exactly,  then  the  angle  A  and  the  time  T  can  be  used 
to  trace  the  signal  trajectory  back  over  its  ray  path 
(called  ray  tracing)  to  determine  the  original  position  of 
the  acoustical  source  [Ref.  1  ]•  However,  full  realization 
of  this  method  in  actual  hydrophonic  tracking  Ls  prevented 
by  two  primary  sources  of  inaccuracies  in  the  process.  The 
first  and  probably  greatest  problem  is  that  tie  speed  of 
sound  profile  can  be  approximated  at  only  a  fe*  locations, 
usually  at  great  cost  to  achieve  even  moderate  accuracy.  In 
addition  the  profile  is  certain  to  flucuate  over  time  and 
location.  The  second  problem,  confounding  the  first,  is 
that  there  may  be  innacuracies,  of  unknown  size,  in  the  time 
data    values   recorded    by   the    sensing    array. 

B.  PURPOSE 

As  noted  in  [Ref.  1]  the  procedure  of  determining  a 
sound  source  position  by  ray  tracing  is  very  sensitive  to 
even  small   errors  in  the  angle   of  elevation   or  speed  of 
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sound  at  the  sensing  array.  Of  course  ray  tracing  is  also 
highly  dependent  on  the  accurate  determination  of  the  the 
transit  time  from  source  to  array.  The  purpose  of  this 
study  is  to  develop  an  appropriate  model  of  the  timing  data 
observed  in  the  hydrophonic  tracking  problem.  Tie  objective 
of  the  desired  model  is  to  estimate  the  error  in  the  time 
values,  and  produce  improved  estimates  of  the  iaitial  angle 
and  ray  transit  time,  so  as  to  reduce  the  ef f 2 ct  of  these 
errors  on  the  ray  tracing  procedure. 

In  pursuit  of  these  objectives  this  thesis  first  reviews 
the  currently  used  models,  which  are  called  the  NAVY  and 
NAVY_A  methods.  Then  four  alternative  models  are  developed, 
called  the  L.S.,  L.S.C.,  M.L.P.  and  M.L.S.  methols.  Finally 
the  performance  of  all  six  models  are  evaluated  and  compared 
using  simulation  studies.  The  comparisons  are  made  under 
the  idealized  conditions  of  a  known  linear  spe=  d  of  sound 
versus  depth  relationship. 

C.   TEACHING  RANGE  CCNFIGOE ATION 

The  tracking  range  which  supplied  data  for  this  study 
consists  of  several  separate  three  dimensional  hydropnone 
arrays  sitting  on  the  sea  bottom.  They  are  lar  ed  out  in  a 
rough  grid,  each  array  being  approximately  750)  feet  from 
the  next,  so  that  the  sound  source  being  tracked  is  never 
more  than  about  5000  feet  from  the  nearest  array.  The 
arrays  are  at  depths  cf  roughly  1000  to  1300  feet. 

Each  array  (see  figure  1.1)  has  four  independent  hydro- 
phones defining  an  orthogonal  coordinate  system.  The  phones 
are  referred  to  as  the  x,y, z,  and  c  hydrophones,  and  are  on 
four  adjacent  vertices  of  a  cube  (see  figure  1.2—  with  sides 
of  length  D  (usually  30  feet).  The  arrays  are  linked  to 
shore  based  computers  by  electronic  cable.  Ths  origin  of 
the  local  coordinate  system  of  each  array  is  at  the  acoustic 
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Figure    1-1         Acoustic   Signal   Detection    by 
Three  Dimensional   Hydrophonic    Array. 


center  of  the  array,  which  is  the  center  of  the  cube  defined 
hy  the  four  hydrophones.  Therefore  the  hydrophones  are  in 
the    positions 


x    : 

y  : 
z    : 

c   : 


(      D    ,  -D   ,    -D  )    /  2 

(  -D    ,  0,-0)72 

(   -D    ,  -D    ,      D    )    /   2 

(    -D,  -0,-0)72 


in   terms   of   the    array's   local   coordinate    system, 
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O* 


(-D,-D,D)/2 


(-D,-D,-Db/2 


/  £-~   "ac°ustic  center 
/^>\  (0,0,0) 


o 


y  :  (-D,D,-D)/2 


x  :  (D,-D,-D)/2 


Figure  1.2    Geometry  of  Hydrophonic  Arnys. 

The  sound  source  is  equipped  with  a  clock  synchronized 
with  the  shore  based  computers  and  emits  a  signil  at  speci- 
fied intervals.  The  times  of  receipt  of  those  signals  at 
the  four  hydrophones  are  recorded,  and  the  corresponding 
travel  times  are  calculated  by  subtraction  of  the  signal 
emission  time. 

D.   APPABENT  POSITIONS 

The  first  step  in  estimating  the  position  of  a  sound 
source  is  to  do  so  under  the  assumption  that  th>  sound  wave 
travelled  its  entire  trajectory  through  water  which  had  a 
constant  speed  of  sound.  The  result,  called  the  apparent 
position,   is  obviously  erroneous,   but  is  then  corrected  by 
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the  ray  tracing  procedure  (a  description  of  wiich  follows 
later)..  The  constant  velocity  value  used  is  usually  the 
estimated  speed  of  sound  V  at  the  depth  of  : he  sensing 
array. 

The  constant  velocity  assumption  implies  that  the  wave- 
front  is  an  expanding  sphere.  Therefore  the  apparent  posi- 
tion is  calculated  using  simple  spherical  equations 
involving  squared  distance  calculations.  Specifically,  if 
Tx  is  the  travel  time  recorded  for  the  x-phone,  then  V«Tx 
is  the  distance  between  the  apparent  position  and  the 
x-phone.  This  distance  is  also  equal  to  the  usui 1  geometric 
distance  between  the  two  positions 

(  X  ,    Y  ,  Z  )        (apparent  position) 
and    (  D  ,-D  ,-D  )  /  2    (x  hydrophone  position) . 
Equating  these   two  squared   distances,   and   equating  their 
counterparts  for  the  other  three  hydrophones,   tie  equations 
(1.1)  are  obtained. 


(  X  -  D/2  )   +  (  Y  +  D/2  )2  +  (  2  +  D/2  ) 


(  X  +  D/2  )   +  (  Y  -  D/2  )2  +  (  Z  +  D/2  )2   =   vV 

y 

2  2 

=   V  T 

z 


V2T2 

x     (1.1) 


(  X  +  D/2  )2  +  (  Y  +  D/2  )2  +  (  Z  -  D/2  ) 


(  X  +  D/2  )   +  (  Y  +  D/2  )2  +  (  Z  +  D/2  )2   =   V2T2 

c 

Assuming  that  the  times  Tx,  Ty,  Tz,  and  T:  are  known, 
(1.1)  is  a  system  of  four  equations  in  three  unknowns  X,  Y, 
and  Z.  This  d verdetermined  system  will,  in  general,  have  no 
exact  solution.  In  fact,  even  if  the  time  values  were 
exactly  correct,  the  system  would  still  have  no  exact  solu- 
tion. This  is  because  the  equations  correspond  to  the 
straight  line  ray  paths  due  to  the  constant  velocity  assump- 
tion, whereas  the  time  values  correspond  to  tie  true  ray 
paths  which  are  not  straight  due  to  the  actual  variation  of 
velocity  of  sound  alcng  the  ray  path.  This  is  a  subtle,  but 
very  important  point. 
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To  throw  out  one  of  the  equations  arbitrarily ,  so  as  to 
reduce  the  system  to  three  equations  in  three  unknowns,  is 
to  throw  away  information  from  one  of  the  hydrophones.  The 
psuedo  solution  currently  utilized  is  to  subtract  the  fourth 
equation  from  each  of  the  first  three,  yielding  a  system  of 
three  equations  in  three  unknowns  wnich  allocs  an  exact 
solution  involving  information  from  all  four  hydrophones. 
However  that  solution  will  not,  in  general,  satisfy  any  of 
the  original  four  equations,  and  is  only  one  of  many  reaorr- 
able    ways   to  choose    an  approximate   solution. 

E.       INITIAL    ELEVATION    ANGLE    AND    RAY    TBANSIT    TIME 

Assuming  that  a  solution  (Xa,Ya,Za)  has  beei  determined 
for  the  apparent  position,  then  the  initial  angle  of  eleva- 
tion is  just  the  angle  of  elevation  of  that  solution,  given 
by     (1.2). 


Al      "      arcsin       (   Za/-\/Xa   +  Ya   +   Za 


The  objective  is  to  find  an  apparent  positioi  (Xa,Ya,Za) 
such  that  (1-2)  computes  an  angle  which  approximates  the 
physically  correct  elevation  angle  as  closely  as  possible. 
The  solution  (Xa,Ya,Za)  and  the  times  Tx,  Ty,  Tz  and  Tc  can 
then  be  used  to  determine  an  appropriate  value  f ) r  T,  which 
is  the  'ray  transit  time',  or  time  of  arrival  of  the  sound 
wave      at    the     acoustic   center      of    the      sensing    array.  The 

currently  employed  method  uses  the  proportional  relationship 
of  equation  (1.3),  where  B  and  Re  are  the  ranj es  from  the 
apparent  position  to  the  acoustic  center  and  to  the  c  hydro- 
phone  respectively,    as  in    (1.4). 
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R 

c 

R 

■  '     ■■    - 

=         V 

T 

1 

T 

c 

/    2 

2         „2 

R      =    - 

v     X      + 

y     +  z 

\/    a 

a          a 

or  T      = 


T      R     \  (1  .3) 

c  *  v 


R 

c 


(1.4) 


R      =      _  /(X    +D/2)2    +     (Y    +D/2)2    +     (Z    +D/2)2 

C  \  /         3  d  3. 


F.       RAY    TRACING 

Whichever  method  is  selected  to  produce  tie  apparent 
position  (Xa,Ya,Za) ,  it  is  transformed  into  the  estimate  of 
the  true  sound  source  position  (X,  Y,Z)  by  the  procedure  of 
ray  tracing.  When  there  is  velocity  layering  ii  the  water, 
the  ray  path  is  no  longer  a  straight  line.  This  is  treated 
using  repeated  applications  of  Sneli's  Law  [Ref.  2  p.  131], 
starting  with  the  layer  of  water  in  which  the  array  sits, 
and  backtracking  upwards  through  successive  velo;  ity  layers, 
until   the   estimated    ray   transit   time   T   is    consumed. 

The  layering  effect  is  artificially  induced  by  the  limi- 
tation that  the  speed  of  sound  can  be  estimated  at  only  a 
finite  number  of  depths,  the  result  of  which  is  commonly 
called  the  water  column.  For  example,  at  the  tracking  range 
studied  the  speed  of  sound  is  measured  every  25  feet, 
starting  at  the  depth  of  12.5  feet.  Hence,  for  example,  the 
third  layer  from  50  to  75  feet  deep  is  assumel  to  have  a 
constant   speed    of    sound   equal    to    that   measured   at    62.5    feet. 

The  first  layer  processed  [Ref.  3  p.  4]  is  the  partial 
layer  lying  between  the  array  and  the  deepest  laf er  boundary 
that  is  shallower  than  the  array,  with  thickness  Z1  (see 
figure    1.3). 

The  incremental  slant  range  in  the  first  Layer  is  S1 
given  by  (1.5),  where  A1  is  fchs  initial  elevation  angle 
estimate. 


16 


sea   surface 
Layer  4 

Ray  Traced 
Position 

Layer  3                          Y 

/v3 

A            "4 

Apparent 
i^-i      Position 

.               yO 
y 
y 

y 

'                 y 

•   y 

/» 

Laver  2                bo/   1  *v 

/  s\      i 

y                                        < 

s 
y 

laver 

Layer  1     v^       I  z 

u 

Acoustic     i 
Center 

boundaries 
sea   bottom 

T» r-  -     /     f  •/  •••/■ f    /     s'    /     /      /'//"/> 

'      .'       '     S          's/'S,''',''',' 

•  s    '  .'  *    '    '    '  '    '   '    '   '    '    '   '  '  ,' 

■' y     '    ,     '    s     s     .     ,     '     '     '    '      '     '    '     '   , 
','    ■    y    <   y    *     -    y   /,''<'  .'  .'  .'  <     ' 

s          '         '        S       *         r          '        ' 
,' ,'        '        '           '       '          '    / 

,    y    '     '  '  s    '    ' 

Figure  1.3    Bay  Tracing  an  Apparent  Position. 


Z1  /  sin  (A1) 


(1.5) 


The  incremental  travel  time  in  the  first  Layer  is  T1 
given  by  (1  .6) ,  where  V1  is  the  velocity  estimated  for  the 
layer    in   which    the    array    is    situated. 
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Tl  =  s!  /  vx  (1.6) 

The  incremental  horizontal  distance  travellel  by  the  ray 
in  the  first  layer  is  H1  given  by  (1.7). 

Hx  =  s1  cos(a1)  (1.7) 

To  determine  the  angle  of  elevation  in  the  next  layer, 
Snell's  law  (1.8)  is  applied,  where  V2  is  the  speed  cf  sound 
estimated   for   that    layer. 


cos (A  )  cos (A  ) 


(1.8) 


Vl  V2 


When  (1.8)    is  solved   for  the  cosine   of  tie   angle  of 
entry  into  the  next  layer,  (1.9)  is  obtained. 

V2  cos (A  ) 
cos(A2)      =      (1.9) 

Vl 
The  procedure  of  computing  the  incremental  values  of 
slant  range,  time  and  horizontal  distance  are  nov  repeated 
for  the  second  layer.  The  overall  procedure  is  repeated 
upwards  through  layers  2,...,n  ,  where  n  is  the  first  layer 
in  which  the  sum  of  the  incremental  travel  times  exceeds  the 
total   ray    transit   time,    as    in    (1.10). 


(Tx    +   T2    +    .     .     .    +    Tn) 


(1.10) 


In  the  last,  uppermost  layer  the  values  Tn,  Sn,  Hn,  and 
Zn  must  be  adjusted  to  compensate  for  overshooting  the  total 
time  T.  The  values  Hi  and  Zi  (i=1,...,n)  are  then  accumu- 
lated  to   get    (1.  11)  . 

Now    the    raytraced   position   estimate    is   given    by     (1.12). 
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a  I         a         a  (1.  12) 


X      =       (    X,    H    )    HxZ    +    Y2 


z  =  z 
The  sensing  array  is  usually  not  aligned  with  the  coor- 
dinate system  of  the  overall  tracking  range.  Therefore  the 
apparent  position,  which  is  in  terms  of  the  local  array 
coordinates,  must  be  changed  by  a  suitable  geometric  trans- 
formation prior  to  ray  tracing  so  as  to  account  for  the 
angle  of  tilt  at  which  the  array  sits  on  the  sea  bottom. 
After  ray  tracing,  the  position  estimate  must  be  again 
transformed  to  account  for  rotation  of  the  array  about  its  Z 
axis  away  from  a  position  which  is  aligned  with  the  range 
coordinate  axes.  Finally  a  simple  translation  must  be 
applied  to  reference  the  position  estimate  to  the  range 
coordinate  system  origin  vice  the  acoustic  ceiter  of  the 
array.  The  end  result  is  a  position  in  terms  of  the  overall 
range      coordinate  system.  These      transformations    are      not 

given  here  because  they  are  used  after  the  estimation  of  the 
initial  angle  and  time,  and  hence  do  not  affect  the  accuracy 
of    those  estimates.       See   [Ref-    3]   for    further   details. 

g.     Discussion 

If  the  velocity  versus  depth  relationship  is  smooth  and 
estimated  accurately,  then  the  ray  tracing  procedure  is 
surprisingly  robust  with  respect  to  the  thickness  of  the 
layers.  For  example,  if  velocity  is  linear  versus  depth, 
and  is  known  exactly,  then  the  exact  hydrophone  times,  ray 
transit  time  and  initial  angles  can  be  computed  [Ref.  4]  for 
any      given    sound      source  position.  Then      the   ray      tracing 

procedure,  with  layers  as  thick  as  25  feet  and  targets  as 
far  away  as  3000  feet,  still  estimates  positions  to  within 
inches  of  each  of  the  true  coordinate  values.  This  seems  to 
indicate  that   errors    resulting      from    position   estimation   are 
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not  due  to  the  approximation  by  layers,  except  pjssibly  when 
there  are  radical  changes  in  the  velocity  pattern  within 
single  layers  such  as  frequently  occur  in  layers  near  the 
water  surface.  Rather  such  errors  apparently  are  due  mostly 
to  inaccuracies  either  in  the  estimation  of  the  speed  of 
sound  profile  itself,  or  in  the  initial  angle  and  transit 
time  estimates.  This  study  will  focus  on  those  arrors  which 
are  involved  in  the  time  values  observed  at  the  hydrophones, 
and  attempt  to  produce  methods  of  initial  angle  and  transit 
time  estimation  which  reduce  the  effects  of  those  errors  on 
the  overall  position  estimation  procedure. 
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II.     COBRENTLY    EMPLOYED    METHODS 

A.       BASIC   CONCEPTS 

The  method  currently  used  for  estimation  of  an  initial 
angle  and  ray  path  transit  time  focuses  on  the  four  spher- 
ical equations  (1.1)  given  in  Chapter  I.  A;  previously 
discussed,    these   equations    have   no    exact    solution    because; 

1.  they  form  an  overdetermined  system  of  four  equations 
in   three    unknowns; 

2.  the  time  values  recorded  at  the  hydrophones  may  be 
innacurate  due  to  unknown  sources  of  error  in  the 
observation    process;    and 

3.  even  if  the  time  values  were  exact,  thef  correspond 
to  a  nonconstant  velocity  profile,  and  so  will  not  be 
be  correct  for  the  constant  velocity  geometry  (spher- 
ical  wavefront)    assumed   by    the   equations. 

As  previously  mentioned,  the  psuedo  solutiDn  chosen  by 
the  current  method  is  to  subtract  the  fourth  spherical  equa- 
tion from  each  of  the  first  three,  and  solve  the  resulting 
system   of   three      equations    in   three    unknowns.  This    method 

has  the  beneficial  quality  that  information  is  rstained  from 
all  four  hydrophones,  whereas  to  just  drop  the  fourth  equa- 
tion (or  any  one  of  the  equations)  without  the  initial 
subtraction  would  cause  complete  loss  of  the  icformation 
from  the  data  recorded  by  one  of  the  hydrophones  .  However 
it  is  important  to  realize  that  tne  solution  thus  obtained 
does  not  actually  satisfy  any  of  the  original  four 
equations. 

It  should  be  noted  that  the  development  of  ti is  solution 
in  [Eef.  3]  is  done  entirely  from  a  geometrical  point  of 
view,       and      does   not      mention    the      system    of      f oi r    spherical 
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equations.  The  text  of  [Ref.  3]  does  not  draw  i  ttention  to 
the  fact  that  the  solution  developed  is  just  one  of  many 
plausible  choices,  none  of  which  will  satisfy  all  four 
spherical   constraints.  Therefore    the      solutioi    chosen     is 

treated  as  though  it  were  the  exact  solution,  only  subject 
to   errors  in  the   observed   hydrophone   time    values.  However, 

even  with  exactly  correct  time  values,  this  currently 
employed  method  will  not  yield  the  true  elevation  angle  and 
ray      transit  time.  This    is      due      to    the      assumption    of      a 

constant  velocity  vice  the  true  nonconstant  velocity 
profile.  This  conflict  introduces  an  automatic  bias  in  the 
initial  angle  and  time  estimates  currently  used  for  ray 
tracing. 

B.       COMPUTATIONS 

To  simplify  notation,  let  (X,Y,Z)  be  the  coDrdinates  of 
the  apparent  position  which  was  formerly  denoted  (Xa,Ya,Za). 
Then  when  the  current  solution  is  applied,  the  first  step  is 
to  subtract  the  fourth  equation  from  eacn  of  the  other 
three,    which  produces  the   equations     (2.1). 

(  X  -  D/2  )2   -   (  X  +  D/2  )2   =   V2  (  T2  -  T2  ) 

c    *       (2.1) 

(  Y  -  D/2  )2   -   (  Y  +  D/2  )2   =   V2  (  T2  -  T2  ) 

c     y 

(  Z  -  D/2  )2   -   (Z+D/2)2   =   V2  (  T2  -  T2  ) 

c     z 

The   solution   to    these  are    easily   oDtained,    as   in    (2.2). 

X      =      V2     (    T2    -    T2    )    /    2    D 

x  (2.2) 

2  2  2  K*-"*-) 

Y      "      V       (    T^    -    T2y    )    /   2    D 

2  2? 

Z      =      V       (    T      -   T      )    /   2    D 

c  z 

Then   the  initial    elevation    angle    estimate    is     (2.3). 


A     =       a 


resin    I   Z  j ~\j    X      +   Y^    +   Z      J  (2.3) 
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The   ray  path   transit  time   to  the   acousti:  center   is 
(2.  4)  ,    where  R  and  Re  are  as  defined  in  Chapter  I  by  (1.4). 

T      =      Tc    R  /   Rc  (2.4) 

This  method  -of  determining  the  apparent  position  shall 
hereafter  be  referred  to  as  the  'Navy  Method',  or  'NAVY'  for 
short . 

C.       ADJUSTMENTS    TO    THE   ORIGINAL    SOLUTION 

Experience  has  shown  that  the  NAVY  method  produces  an 
apparent  position  estimate  which  usually  can  be  improved  by 
an    adjustment    which    is   described    in    this    section,. 

The  cosine  of  the  angle  between  the  i-th  axis  and  the 
straight  line  from  the  origin  out  to  the  apparent  position 
is  called  the  i-th  direction  cosine  Ci.  It  is  a  fact  of 
geometry  that  the  sum  of  the  squares  of  the  three  direction 
cosines  must  equal  unity.  Therefore  the  method  is  now 
adjusted  to  reflect    that  constraint. 

The  direction  cosines  used  are  the, angles  nade  by  the 
ray  path  at  the  c  hydrophone.  Therefore  the  (X,Y,Z)  coordi- 
nates calculated  by  the  original  method  are  first  translated 
to  coordinates  referenced  to  the  c-phone  as  the  temporary 
origin,    as   in    (2.5). 

X      =   X   +   D/2  Y      =   Y   +   D/2  Z      =    Z   +   D/2  _     ,. 

c  '  c  ■  c  (2.b; 

Therefore   the   three   direction  cosines  ar=   given  by 
(2.6)  . 

Cv      -      X      /   V  T  C         =      Y      /   V   T  C         =      Z      /   V   T  ,0     ,. 

x  c  c  y  c  c  z  c  c(2.6) 

The  denominators  in  (2.6)  are  all  V«Tc  beci  use  that  is 
the  range  from  the  apparent  position  to  the  c  hydrophone,  as 
estimated   by   the    time  from    the   c    hydrophone.         Ideally    these 
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cosines    should    add    to   unity    when    squared.         Tnersfore    if   DCC 
is      the    'direction      cosines      correction*       factor    defined      by 
(2.7)  ,    then   DCC    should   be   close    to    one. 


^  /  C2  +  c2  +  c2 
x  y  z 


(2.7) 


Deviation  of  DCC  from  unity  is  interpreted  as  an  indica- 
tion of  receiver  timing  errors,  array  malformation  or 
invalid  data  at  one  or  more  of  the  hydrophones  [Ref.  3 
p.C-3].  Currently      if      DCC        lies      outside      tie      interval 

(0.98,1.02),  the  data  is  discarded  as  being  excessively  full 
of  error.  The  direction  cosines  of  the  remaining  acceptable 
data  points  are  rescaled  using  (2.8)  to  assure  satisfaction 
of   the   direction   cosines  constraint. 


Cx     =     Cx  /  DCC  Cy     =     cy  /  DCC  cz     =     c2  /  DCC         (2.8) 

A      corrected   set      of      new      coordinates   are      ;  omputed     by 
(2.9),    still  being    referenced    to    the   c   nydrophone. 

X      =      V  T      C  Y      =      VTC  Z=VTC 

c  cx  c  cy  c  cz  (2.9) 

These  are  then  translated  by  (2.10)  to  coordinates 
referenced  to  the  acoustic  center. 

X    =    Xc-D/2  Y    =   Yc-D/2  Z    =    ^    -    D/2  ^^ 

This  adjusted  method  of  determining  the  apparent  posi- 
tion shall  hereafter  be  referred  to  as  the  'Navy  Adjusted 
Method',    or    'NAVY_A'    for   short. 


D.       DISCUSSION 

When  a  sound  source  is  within  the  detection  range  of 
more  than  one  sensing  array,  each  array  produces  timing 
data.         The   data    from  each    array    can    be    processed    to   produce 
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a  position  estimate.  Ideally  these  multiple  estimates  of 
position  will  be  in  reasonably  close  agreement.  However 
experience  with  actual  tracking  data  has  shown  that  this  is 
not  the  case  in  many  of  the  multiple  detection  opportuni- 
ties. This  tendency  toward  disagreement  between  multiple 
estimates  of  the  same  position  is  commonly  called  the  cross- 
over, or  crosstalk,  problem.  This  problem  often  occurs  when 
the  sound  source  is  moving  away  from  the  tracking  domain  of 
one  array  into  the  tracking  domain  of  another.  This  study 
focuses  on  improvement  of  the  initial  angle  and  time  esti- 
mates, which  hopefully  will  help  alleviate  the  crossover 
problem. 

The  current  choice  of  a  'best*  compromise  solution 
appears  to  be  based  on  reasons  of  simplifying  geometry  and 
calculations.  These  are  worthwhile  objectives,  but  do  not 
in  themselves  reflect  the  need  to  estimate  accurately  the 
initial  elevation  angle  and  ray  transit  time.  Since  there 
exist  physically  correct  values  for  both  the  ai gle  and  the 
time,  those  values  will  produce  the  exact  position  after  ray 
tracing,  provided  that  the  velocity  -profile  is  known 
exactly.  The  desire  then  is  to  estimate  these  true  values 
as   accurately   as  possible. 

The  question  at  this  point  is  whether  or  not  the  direc- 
tion cosines  adjustment  causes  the  estimated  apparent  posi- 
tion to  fce  closer  to  the  true  apparent  position.  Experience 
has  indicated  that  it  does  [Bef.  3  p.C-7].  However  the 
effect  of  the  adjustment  can  be  interpreted  in  terms  of  the 
original  four  spherical  equations  (1.1).  The  rescaling  of 
the  direction  cosines  given  by  (2.6),  so  as  to  assure  that 
their  squares  add  to  unity,  is  equivalent  to  rescaling  the 
quantities  in  (2.5)  so  as  to  assure  that  their  squares  add 
to  (V«Tc)2.  That  is  exactly  the  constraint  stated  by  the 
fourth  spherical  equation  of  (1.1).  So  the  effect  of  the 
adjustment      is        to      require        that      the        fourti       equation, 
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concerning  the  data  at  hydrophone  c,  is  exactly  satisfied. 
This  requirement  will,  in  general,  assure  that  the  other 
three  equations  are  not  satisfied.  Since  experience  shows 
that  the  adjustment  often  improves  the  solution,  this  seems 
to  imply  that  the  fourth  eguation  is  somehow  more  important 
than  the  other  three.  Or  it  may  just  be  that  exact  satis- 
faction of  one  of  the  equations  usually  assures  a  reasonably 
good    compromise    solution. 

Tc  summarize,  the  NAVY  method  provides  a  useable 
apparent  position  suitable  as  input  for  ray  tricing.  But 
the  direction  cosines  adjustment  used  in  the  NAVY_A  method, 
for  reasons  not  understood  at  this  time,  appears  to  improve 
that  position  as  indicated  by  test  results.  Ta ose  results 
are  supported  by  the  results  of  this  thesis  (see  Chapter  V)  . 
However,  as  will  be  demonstrated  by  the  example  considered 
in  the  next  section,  the  DZC  correction  factor  o:  the  NAVY_A 
method  has  an  effect  which  must  be  something  more  than  just 
the    smoothing   of   timing    errors. 

E.       A    COMPUTATIONAL    EXAMPLE 

Eor  the  purposes  of  illustration  and  comparison,  suppose 
that  a  30  foot  sensing  array  is  at  a  depth  of  1300  feet, 
that  the  coordinate  system  origin  is  at  the  array  acoustic 
center,  and  that  the  array  arms  are  parallel  to  the  coordi- 
nate axes.  This  implies  that  the  four  hydrophones  are  in 
the    positions: 

x    :        (    15    ,    -15    ,    -15) 

y    :        (-15    ,       15    ,    -15) 

z    :        (-15    ,    -15    ,       15) 

c    :       (-15    ,    -15    ,    -15) 
If   there  is    a   sound   source   known    to   be   in   position 

(    1000    ,    3000    ,    900    ) 
then    the   depth    of   that   source    is    1300    -    900    =    40)       feet. 
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Figure  2.  1    Sample  Depth  Versus  Velocity  Profile. 
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Figure  2.1  shows  the  estimated  sound  velocity  profile 
which  was  estimated  for  the  data  used  in  the  course  of  this 
thesis.  As  :an  be  seen,  the  profile  is  primarily  linear  at 
depths  greater  than  100  feet.  The  profile  at  depths  greater 
than  100  feet  is  reasonably  approximated  by  the  linear 
relationship 

V      =      4840.7      +       0.03314    •    DEPTH    . 

Therefore  suppose  that  in  this  example  problem  the 
velocity  profile  is  known  exactly,  and  is  given  Dy  the  above 
linear   relationship. 

Under  these  circumstances,  witn  known  lins ar  velocity 
profile,  and  known  sound  source  location,  the  exact  times  of 
arrival  of  the  sound  wave  at  the  four  hydrophones  can  be 
computed  using  the  methods  set  forth  in  Appendix  A.  Those 
exact    times    (in    seconds)    are: 

Tx    :       0.6779686893  Ty    :       0.674  2257788 

Tz    :       0.6782243197  Tc    :       0.6798324156    . 

The  corresponding  exact  values  for  the  initial  elevation 
angle,  ray  transit  time  and  resulting  apparent  position  are 
also  directly  computable.  Those  computations  will  hereafter 
be  be  known  as  the  EXACT  method,  and  will  produce  the 
correct  true  position  after  ray  tra-cing.  The  results  of  the 
EXACT  method  are  given  in  Table  I,  along  with  the  corre- 
sponding apparent  position  estimates  produced  waen  the  two 
methods,  NAVY  and  NAVY_A,  are  applied  to  the  (errorless) 
time    values. 

At  first  glance  the  differences  in  Table  I  might  seem 
rather  small.  However  it  is  important  to  recall  that  these 
are  produced  under  the  ideal  conditions  of  a  very  smooth  and 
exactly  known  velocity  profile.  These  idealizations  are  far 
from  the  realities  of  a  nonlinear  velocity  profile  which  is 
estimated  by  a  procedure  involving  errors  which  are  unknown 
and  probably  significant.  Such  realities  mighc  well  cause 
the    differences    in    Table  I    to    be   significantly   larger.         The 


28 


TABLE    I 

Single   Exam 
NAVY   and 

pie  Comparison   of 
HAVY_A    Hethods 

Method 

Transit    Time 
T    Tsecs) 

Slev.    Angle 

NAVY 

0.67526969 

15. 26459 

NAVY_A 

0. 67527027 

15.26461 

EXACT 

0. 67527043 

15.27002 

Apparent 

£2§i.tion        Estimate 

NAVY 

(    1005.957 

/ 

3017.874    ,       868.143    ) 

NAVY_A 

(    1006.087 

r 

3018.259    ,       868.255    ) 

EXACT 

(    1005.966 

* 

3017.899    ,       868.555    ) 

_ 

nature  and  size  of  those  differences  remain  difficult  to 
determine  until  more  is  known  about  the  the  velocity  profile 
estimation  errors  and  their  effect  on  the  position  esti- 
mating process.  In  any  event  even  the  small  differences  in 
Table  I  might  be  magnified  during  the  ray  tracing  process 
under   the   conditions   cf   a  realistic   velocity   profile. 

The  differences  illustrate  the  very  important  point  that 
the  direction  cosines  adjustment  causes  changes  L_n  the  esti- 
mates even  when  the  time  data  is  free  of  all  er:_or.  Hence 
the  deviation  from  1.0  of  the  correction  factor  DCC  is  not 
just  due  to  array  malfunction,  receiver  timii g  error  or 
other   sources   of   non-valid    data   as   previously   assumed. 

In  the  example  above,  the  NAVY_A  method  produced  a 
slightly  better  time  and  angle  value  than  the  i  AVY  method. 
However,  in  this  same  example  the  NAVY_A  method  produced  an 
apparent  position  estimate  which  is  slightly  farther  away 
from    the      EXACT    answer      tnan    the      position    estimated      by    the 
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NAVY  method.  This  is  only  one  example,  and  its  results 
should  not  be  generalized.  However  it  illustrate  s  the  point 
that  apparently  the  true  effect  of  the  adjustment  may  not  be 
well  understood. 

Furthermore,  the  adjustment  seems  to  place  heavy 
emphasis  on  the  time  value  recorded  at  hydrophone  c. 
Therefore  the  effect  of  the  adjustment  may  well  depend 
largely  on  the  accuracy  of  that  one  particular  data  value, 
which  is  a  relatively  unbalanced  dependence  in  the  presence 
of  data  errors. 
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III.  PLANAR  BA7EFR0NT  MODELS 

A.  TEE  PLANAR  WAVEFRONT  ASSUMPTION 

When  a  sound  wave  travels  in  constant  veloci:  y  water,  it 
expands  in  the  shape  of  a  sphere.  If  the  velocity  profile 
is  variable  instead,  but  is  reasonably  well  behaved  versus 
depth,  then  the  expanding  wave  is  a  smooth  distortion  of  a 
spherical  surface.  In  either  case,  if  the  wave  has 
travelled  a  long  distance  when  it  arrives  at  a  hydrophonic 
array,  then  that  small  piece  of  the  wavefront  # hich  passes 
through  the  30  foot  cube  spanned  by  the  array  may  be  approx- 
imated reasonably  by  a  flat  planar  surface.  This  approxima- 
tion is  the  basis  of  the  planar  wavefront  models  developed 
in  this  chapter. 

B.  PLANE    EQUATIONS 

A  plane  in  space  is  fully  defined  by  idei  tif ying  any 
point  (X0,Y0,Z0)  on  the  plane,,  and  also  a  vector  (C1,C2,C3) 
of  unit  length  which  is  perpendicular  to  that  plane.  The 
vector  is  called  the  unit  normal  vector  for  that  plane. 
Then  any  point  (X, Z,Z)  on  that  plane  must  satisfy  the  equa- 
tion   of   the   plane,    namely    (3.1). 

cx   (  x  -  xQ  )   +  c2    (  y  -  yq  )   +  c3    (  z  -  zQ  )     =     o      (3-1J 

The  perpendicular  distance   from  the  plane  t)   any  point 
(X1,Y1#Z1)  not  on  the  plane  is  the  absolute  value  of  (3.2). 

Cl  (  Xl  "  X0  >  +  C2  (  Yl  "  Y0  }  +  C3  (  21  "  Z0  >      (3.2) 
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C-   BASIC  HODEL  EQUATIONS 

Consider  a  coordinate  system  whose  origin  is  at  hydro- 
phone c,  and  whose  axes  are  aligned  with  the  three  array 
arms.  Let  C1,  C2  and  C3  be  the  direction  cosines  on  the  X, 
Y  and  Z  axes  respectively  for  the  vector  from  the  origin  to 
the  apparent  sound  source  position.  Then  (C1,C2,C3)  is 
itself  a  vector,  of  unit  length,  which  is  perpendicular  to 
the  planar  wavefront  emanating  from  the  soi  nd  source. 
Therefore  (C1,C2,C3)  can  be  used  as  the  normal  vector  for 
the  wavefront  plane. 

In  the  coordinate  system  referenced  to  the  c-phone  as 
the  origin,  the  acoustic  center  has  coordinates  (D,D,  D)/2, 
where  D  is  the  length  of  an  array  arm.  When  the  soundwave 
plane  arrives  at  the  acoustic  center,  it  will  have  the 
equation  (3  .3)  . 

CX      +      C      Y      +      C      Z      =      D    (    C.    +    C_    +    C-    )    /    2 
1  2  3  12  3  (3.3) 

The  x-phone  has  coordinates  (0,0,0),  and  : he  distance 
between  it  and  the  soundwave  plane  at  the  acoustic  center  is 
(3.  4),    which   then   simplifies    to    (3.5). 


Cx    (   D/2    -  D   )    +  C      (    D/2    -   0   )    +  C      (   D/2   -   0   ) 


(3.4) 


D    (    -Cx    +    C2    +    C3    )    /    2  {3m5) 

This  distance  is  measured  in  a  direction  perpendicular 
to  the  wavefront  plane,  and  so  is  measured  in  the  direction 
of    travel   of  the   soundwave.  Therefore    the    same    distance   is 

also    equal    to    (3.6)  . 

V(T*-T)  (3.6, 
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In     (3.6)    V    is   the  velocity   of    sound    at    the   array,      Tx    is 
the    time    of   arrival    of   the    sound    wave    at    the    x-phone,       and  T 
is    the      time   of    arrival      of    the      sound    wave   at      the    acoustic 
center.        The   term    "distance"    is      used   loosely   in     (3.5)       and 
(3.6)  ,      because    these   quantities    may    be    negative.  The   true 

distances  are  the  absolute  values  of  these  quantities. 
Since  the  next  step  is  to  equate  these  two  distal ces,  it  is 
only  necessary  to  show  that  these  two  quantities  always  have 
the  same  sign.  There  are  two  cases  to  consider,  depending 
on  whether  the  first  component  of  the  apparent  position  is 
positive  (X>0),  or  negative  (X<0)  .  Let  (X',0,0)  be  the 
intersection  of  the  X  axis  with  the  wavefront  plane  as  it 
passes  through  the  acoustic  center.  Then  (3.3)  can  be  used 
to   solve   for  X',    namely 

X*      =      D     (C1    +   C2    +    C3)     /    2    C1. 
Now    consider   the    case  where   X>0.      Then    C1>0   ilso,    and   if 
(3.6)     is   positive,    then 

Tx    >    T 
which    implies   that      the   wave  arrives   at      the   acoustic   center 
before   it  arrives  at    the  the   x   hydrophone.      Therefore,    since 
X>0 ,       the   plane    at    the  acoustic      center    will   intersect    the   X 
axis    at    a   pcint    beyond    the    x   hydrophone,    or    X'>D,     and    hence 
X»    >    D         =>       (C1    +    C2    +    C3)/    2    C1         >       1 

=  >      C1    +   C2    +    C3      >    2    C1  (sines     C1>0) 

=  >       -C  1    +    C2    +   C3      >       0 

=  >   (3.5)  is  positive  . 
A  parallel  argument   applies  for  the  case   of  X>0,   thus 
establishing  that  (3.5)  and  (3.6)  always  have  the  same  sign. 
Equation   (3.7)    is   the   result    of   equating   these   two 
quantities. 

-C±    +    C2  +  C3  +  (  2  V  T  /  D  )  -  (  2  V  ?      /  D  )    =   0      (3.7) 
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For  convenience,  let  K  =  2V/D,  and  then  apply  the  same 
logic  to  the  distances  of  the  y,  z,  and  c  hydrophones  from 
the  wavefront  plane  as  it  passes  through  the  acoustic 
center,  to  obtain  the  equations  of  (3.8). 


(3.8) 
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The  system  (3.8)  is  four  equations  in  four  ui knowns,  and 
is  the  planar  model's  version  of  the  equations  (1.1).  The 
unknowns  in  (3.8)  are  C1,  C2,  C3  and  T.  However  there  is 
the  additional  constraint  that  C1,  C2,  and  C3  a:  e  direction 
cosines,  and  therefore  the  direction  cosines  constraint 
(3.9)  is  a  fifth  equation,  creating  a  a  system  of  five  equa- 
tions  in   four   unknowns. 

^2  2  2 

Cl      +      S      +      C3  "      1  (3.9) 

Generally  there  is  no  set  of  values  (C1,C2,:3,T)  which 
will  satisfy  ail  five  equations  at  once.  This  is  because  of 
the  realities  of  a  nonplanar  wavefront  and  the  presence  of 
timing  errors.  The  next  section  developes  a  method  that 
produces  set  of  values  for  the  unknowns  which  is  intended  to 
satisfy  those  equations  reasonably  well. 

D.   HIHIBIZATION  OF  SUM  OF  SQUARED  ERRORS 

Let  Ei#  (i=1,2,3,4)  be  the  value  of  the  left  hand  side 
of  the  i-th  equation  of  (3.8).  Then  Ei  measures  the  amount 
of  error  in  the  i-th  equation  caused  by  the  chosen  solution. 
It  is  impossible  to  have  Ei=0  for  all  i=1,2,3,4.  However 
some  compromise  may  be  made.  Specifically  tha  compromise 
chosen  here  is  the  classic  minimization  of  (3.10) ,  the  sum 
of  squared  errors. 
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^2  „2  2 


(3.  10) 


Minimization  is  to  be  done  subject  to  tb.2  direction 
cosine  constraint  (3-9).  Application  of  the  Lagrange  multi- 
plier technique  [Ref.  5  p. 55]  calls  for  the  minimization  of 
(3.11)  over  all  possible  choices  of  CI,  C2,  C3,  T  and 
lambda. 


k<  -*&<  -\ 


(3.11) 


Taking    the       partial    derivative      of   L      with    rs  spect      to    T 
yields    (3. 12)  . 
3l 


h  t 


Y      (-2K  E.) 
iti 


-2K    (    E±    +   E2    +   E3    +    E4    ) 


(3.  12) 


2    K     [-2K   T      +      K    (    T      +    T      +    T      -    T4    )  j 


Eguating  (3.12)  tc  zero  and  solving  for  T  yields  immedi- 
ately the  appropriate  estimate  (3.13)  of  T,  the  ray  transit 
time. 


(    T      +   T      +   T      -    T       ) 

2  1  2  3  4' 


(3.13) 


The      partial      derivative    of      L      with      respect    to      C1      is 
(3.14)  . 

2     (    E      -    E      -   E      +   E      ).    -  2  \  C 

J  -1  (3.  14) 

4C      +    2KT   +   K    (    T    -T    -T    -T      )       -   ^  Ci 


Be, 


If    (3.13)    is    used  for   T    in    (3.14),    then    (3.15)     results. 


3l 


2  (    4    -   \)    C1        +      2    K    (    Tx    -    T4    )]    (3.  15) 
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If  the  same  procedure  is  used  for  the  partial  deriva- 
tives of  L  with  respect  to  each  of  C1,  C2  and  C3,  and  all 
are    equated   to    zero,    then  the    results    are    (3.  16)  . 


(   4   -\)    C. 


2    K    (    T„    -    T.     ) 
4  l 


i  =    1,2,3 


(3. 16) 


In   order     to   solve  for      the   Lagrange      muitipL ier    lambda, 

square   both   sides  of    the   three   equations   of    (3.16),      and   add 

•the      resulting    equations      together.         Then      use      the    sum      of 

squares      constraint     (3.9)         and   solve      for      iambi  a    to      yield 

(3.17). 


2    K 


■    3 
I    (T. 

.1  =  1 


-    T.) 


(3.17) 


Substitute  (3.17)   in  (3.16)   and  simplify  to  obtain  the 
appropriate  estimate  of  Ci,  namely  (3.18). 


C. 

i 


T„  ~  T. 
4     l 


y  (  t.  -  t.  ) 


i  =  1,2,3 


(3.18J 


The  choice  of  sign  in  (3.18)  is  determined  by  the  fact 
that  Ci  is  positive  if  and  only  if  the  sound  wav*  arrives  at 
the  i-th  phone  before  it  arrives  at  the  c-phone,  which  in 
turn  implies  that  (T4-Ti)>0. 

E.   THE  IEAST  SQUARES  METHOD 

In  summary,  the  first  alternative  method  for  determining 
an  apparent  position  starts  with  estimation  )f  the  ray 
transit  time  to  the  acoustic  center  by  (3.19).  Then  the 
apparent  position  estimates  are  computed  using  (3.20). 


(  T   +  T2  +  T3  -  T4  )  /  2 


(3.19) 
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This  method  shall  be  hereafter  referred  to  as  the  'Least 
Squares  Method1,  or  •L.S. '  for  short.  The  apparent  advan- 
tages   of   the  L.S.    method   are   that: 

1.  all  four  hydrophone  times  have  equal  weight  in  a 
simple  expression  for  the  ray  transit  time  T,  rather 
than  using  an  expression  so  heavily  dependent  on  Tc 
as   in  the   NAVY  and   N AVY_A    methods; 

2.  the  differences  of  squared  time  values  t  hich  appear 
in  the  solutions  of  the  NAVY  methods  are  avoided  in 
the  L.S.  method,  thereby  lessening  tie  tendency 
toward  computational    roundoff    problems; 

3.  the  direction  cosines  already  add  to  unity  when 
squared,    requiring    no   arbitrary  adjustmen:  ;    and 

4.  computation  of  the  initial  angle  allows  cancellation 
of  several   terms,      resulting   in    the    simple    expression 

(3.21)  . 


arcsm 


(T4    "T3    '/VZ/^    "Tj    ''J       (3-21) 


F.       BIAS    IN    THE    LEAST   SQUARES    METHOD 

Unfortunately  a  potentially  serious  problem  exists  with 
the  L.S.  method.  That  concerns  the  consequeices  of  the 
assumption   of   a    planar   wavefront.  The    effect    is    difficult 
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to  derive  explicitly  because  it  is  difficult  : o  determine 
the  true  shape  of  a  wavefront  corresponding  to  a  nonconstant 
velocity  profile.  However  if  it  is  assumed  that  a  spherical 
assumption  is  more  accurate  the  planar  assumption,  then  the 
bias  can  be  estimated  roughly,  and  then  subtracted  from  the 
original  L.S.  position  estimate.  This  is  still  a  difficult 
problem  because,  as  previously  discussed,  the  four  basic 
spherical   equations    themselves   have   no  exact   solution. 

Nevertheless,  as  a  rough  estimate  of  the  bias,  the 
following  procedure  is  used.  First  estimate  :  he  apparent 
position   by   the    L.S.      method.  Then   calculate    the   straight 

line  distances  from  that  position  to  each  of  tha  four  array 
hydrophones.  Divide   those      distances  by      the    velocity      of 

sound  at  the  array  to  obtain  the  corresponding  times.  Use 
these  times  to  recalculate  the  apparent  position  using  the 
L.S.  method  again.  The  difference  between  the  original  and 
recalcualted  L.S.  positions  roughly  measures  the  error  that 
would  re  made  by  the  L.S.  method  when  it  is  applied  to  the 
time  values  which  correspond  to  a  spherically  spreading 
sound  wave  whose  source  is  in  the  vicinity  of  the  original 
apparent  position.  Therefore  this  difference  caa  be  used  as 
a  bias  vector  which  can  be  subtracted  from  the  original  L.S. 
solution.  This  bias  correction  is  the  basis  of  the  method 
set    forth  in  the   next  section. 

G.       THE    LEAST   SQUARES   CORRECTED    METHOD 

The  second  alternative  method  for  estimating  an  apparent 
position   is   as    follows: 

1.  calculate  the  apparent  position  P1  by  using  the  L.S. 
method; 

2.  calculate  the  distances  from  that  position  to  the 
four  hydrophones,  and  convert  them  to  times  by 
dividing    by    the   speed    of   sound    at    the    array; 
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3. 


a. 


5. 


use  the  new  times  to  calculate  a  new  apparent  posi- 
tion  P2,    using  the   L.S.    method; 

calculate  the  difference  vector  P2-P1  (see  figure 
3.1),  and  sultract  it  from  the  original  position  P1 
to   obtain    the    corrected    position   P 

P      =      P1    -    (    P2   -    P1    )       =      2   PI    -    P2 
finally    adjust   the   raj    transit      time   T    calculated    for 
the    original    position -P1,         by    using   the    proportional 
transformation: 

T»    =    T    *    R    /    EM 
where  R    is   the  range    to   the   new   position   ?  ,    and    R1    is 
the    range   to    the   original    position    P1. 


Ji 


Recalculated 

L.S.  Estimate 

"  =  P2  -  PI 


Original  L.S, 

\  Estimate 
\ 


Adjusted  L.S.C 
Estimate 


P  =  PI  -  E 


Array 
Center 


Figure  3.1    Bias  Adjustment  for  the  L.S.  Method. 

This  method  shall  hereafter  be  referred  to  as  the  'Least 
Squares  Corrected  Method1,  or  L.S.C.  for  short.  The  proper- 
ties of  the  L.S.C.  method,  like  those  of  the  Ukf Y_A  method, 
are  not  well  understood  at  this  time.   It  is  offered  only  as 
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an  alternative  which  may  combine  the  beneficial  properties 
of    the   l.S.    and    NAVY    methods,    namely    that    it    will: 

1.  estimate  a  ray  transit  time  value  equally  dependent 
on  all  four  data  times,  thereby  smoothing  out  exces- 
sive error  in    any   one  of    the    time    values;    and 

2.  reflect  the  spherical  wave  assumption,  believed  to 
provide  a  more  accurate  description  than  the  planar 
assumption. 

For  targets  at  a  range  of  3000  feet,  the  length  of  the 
error  vector  P2-P1  varies  from  0  to  as  much  as  10  or  11  feet 
(see  Table  II).  The  error  vector  lengths  seem  to  be  depen- 
dent en  both  azimuth  and  elevation  angles  of  the  target  from 
the  array.  These  patterns  indicate  a  potential  for  further 
investigation   to   relate    estimation   errors    to   suci     variables. 

H.       HAXIHUfl   LIKELIHOOD   CALCULATIONS 

One  shortcoming  of  all  methods  described  thus  far  is 
that  none  of  them  can  be  used  to  produce  an  estimate  of  the 
underlying    error  in      the   time    data   values.  The    method   set 

forth  in  this  section  is  a  first  attempt  to  estimate  that 
noise,  and  is  again  based  on  the  assumption  of  a  planar 
wavef ront. 

Let  Ti  be  the  time  recorded  from  by  i-th  hydrophone. 
Let  Ui  be  the  true  time  which,  under  absolutely  error  free 
conditions,  would  have  been  recorded  at  the  i-th  hydrophone. 
An   assumption   of   Gaussian   noise   is   now    made,    namely    that 

Ti      =       Ui    +    Ei 
where   the   Ei  are      independent   identically    distributed   normal 
random   variables    with  mean    zero   and   variance    32.         Therefore 
the    Ti    (i=1 ,2,3,4)  are   also   normally   distributed      with    the 

same    variance,    but    with   means    Oi     (i  =  1,2,3,4). 
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Now  let  Cj  be  the  j-th  direction  cosine  (j=1,2,3).  Then 
geometrically,  using  the  assumption  of  a  planar  wavefront, 
Cj  is  given  by  (3.22),  where  V  is  the  speed  of  sound  at  the 
array,    and    D  is    the    array   dimension    (30    feet) . 


c. 

3 


V     (    U      -    U.     )    /    D 
4  3 


(3.22) 


letting    U   =    0*4 ,    .then    (3.22)    can    be    solved    for    each    Uj    in 
terms   of   (J   and    Cj,    as  in    (3.23). 


u 


U    -    D    C.    /   V 

3 


(3.23) 


The    probability    density      of   each    time    value      Ti    is    given 
by    (3.24). 

(3.24) 


1 


W   - 


V^f 


exo 


-     (    T.     -    U. 

l  i 


2    S 


Using    (3.23)    in     (3.24),       and   multiplying    the    four   densi- 
ties   together,    the    likelihood    function    (3.25)    is    formed. 


exp 


-1  4 

y.  V   (T.     -    U    +   DC./V) 

1=1 


2    S 


(3.25) 


In  (3.25)  C4  has  been  used  for  notational  ;  onvenience, 
and  is  defined  to  be  zero.  Then  the  log-likelihood  function 
is   formed  by  taking    natural    logs   of    (3.25),    yielding    (3.26). 

1  4  2 

L        =      -2    In  (2  7f  )    -    4    ln(S) V    (T.-U+DC./V)  (3.26) 

1  2    S2      &1      L  i  *  ' 

Since  the  values  of  the  Ci  are  to  be  direction  cosines, 
the  usual  direction  cosines  constraint  must  be  added  to  L1 
to   form    the   Lagrangian  function    L2    given    in    (3.27) . 
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L2      =      Ll    -   X        I    «£       -    1 

-1  =  1 


(3.27) 


The  objective  is  to  maximize  L2  over  tie  possible 
choices  cf  C1,  C2,  C3,  U,  S  and  lambda.  Ignoring  S  for  the 
moment,  maximization  of  L2  can  be  acheived  by  minimization 
of   1    given    in    (3.28). 


L      =    2    S2 


]T    (T.    -    U    +    DC./V)  +^ 


r    3 


1=1 


E  <c  >  :  > 

1=1 


(3. 28) 


Take    partial    derivatives   of    L   to    get    (3.29)    i  nd    (3.30) 
B  L  2D 


3c. 

i 


v     (T.    -    U    +   DC./V)       -      2   \C. 


(3.29) 


Bl 

3u 


r    4 


£ 


-    4U      + 


3 


(3.30) 


Equate  (3.30)  tc  zero  and  solve  for  U  to  obtain  (3.31), 
the  maximum  likelihood  estimate  for  the  time  at  the  c  hydro- 
phone. 


)     T         + 

k  j 


h)/ 


(3.31) 


Equate      the    three      equations    of       (3.29)         to    zero,         and 
multiply  each   equation  by   Ci   respectively,    to    ob:  ain    (3.32). 

\     2  D     /  D     2  \ 

XCi        =     —  (TxCi     "     UCi      +-Ci)  i=   1^2,3      (3.32) 

Add   the   three   equations    of     (3.32)     together,       and    use    the 
direction   cosine    constraint    to    obtain     (3.33). 

r      3  3  -j 

X    =    —       X>iCi     "     U     ICi     +-V"  (3.33) 

_i=l  i=l 

Equation    (3.34)     results    when     (3.29)     is    equated    to   zero. 
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-    u 


c. 

1 


i   =   1,2,3 


(3.34) 


£X- 


D   '*•         V 

Substitute  (3.33)  into  (3.34)  to  outain  (3.35),  the 
maximum  likelihood  estimate  of  the  direction  cosines,  where 
U  is    the   estimate  calculated  by    (3.31). 

u 
(3.35) 


T.       - 

l 


3 

y  t.c. 


-    u    V  c. 


£i  3 


As  the  reader  will  perhaps   have  noticed,   tie  equations 


of  (3.35)    define  each  of  the   unknowns  Ci  in  terms 


all 


three  unknowns.  Such  a  structure  suggests  that  (3.35)  can 
be  used  as  an  iteration  function.  That  is,  i:  reasonable 
initial  values  are  used  for  the  three  unknowns  in  the  right 
hand  side  of  (3.35)  then  new  values  are  produced.  Repeat 
the  process  using  the  new  values  until  the  answer  stabilizes 
within  acceptable  tolerances.  Although  convergence  to  the 
correct  solution  is  not  guaranteed,  the  method  has  never 
failed  for  the  equations  of  (3.35).  Unlike  the  L.S.  method, 
(3.35)    does  not    have   any   known  closed    form   solution. 

Returning    to    the    standard    deviation    S,      take    the    partial 
derivative    of    (3.27)     with   respect    to   S    to    get    (3.  36) . 


ds  s  s3     ki  x  v     1 

Multiply   (3.36)   by   S3  and  solve   for  32 

maximum   likelihood   estimate   of  the   variance, 
(3.37)  . 

i(Ti-u)  +-fici 


1=1 


i=l 


(3.36) 

to   get      the 
given      by 


(3.37) 
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I.       THE    HAXIMaa    LIKEIIHOOD    PLANAfi    METHOD 

In      summary,       the      third      method      for    estimation      of      an 
apparent   position   is    as   follows: 

1.  let    0  =    T4    initially; 

2.  use  the  L.S.  method  to  calculate  the  initial  values 
for    Ci,     (i=1,2,3),    using    (Se321); 

3.  use  0  and  Ci  (1=1,2,3)  in  the  right  hand  side  of 
(3.35)    to   obtain   new    estimates    for    Ci    (i  =  1  ,2,3)  ; 

4.  recalculate    U,    using     (3.31); 

5.  reiterate  steps  3  and  4  until  tne  values  Ci  (i=1,2,3) 
and    U  converge   within  acceptable    tolerancss; 

6.  calculate    S2,    using     (3.37)  ; 

7.  calculate  the  estimated  apparent  position  relative  to 
the    c-phone,     using    (3.38); 

Xc      =      VUC1  Yc      =      VU    C2  Zc      =      VUC3     (3>38) 

8.  lastly  translate  this  solution  and  its  o  rresponding 
time  estimate  to  a  solution  and  time  relative  to  the 
acoustic  center,  using  (3.39),  where  E  and  Re  are  as 
defined  by  (1-4)  in  Chapter  I. 


X      =      X      +    D/2  Y      =      Y      +   D/2  Z    =    Z      +    D/2 

c  c  c 

T      =      U    R  /    R 
c 


(3.39) 


This  method  shall  be  hereafter  referred  to  as  the 
'Maximum  likelihood  Planar  Method',  or  M.l.P.  for  short. 
Originally  the  hopes  for  this  method  were  rather  high,  espe- 
cially since  it  was  the  first  method  to  produce  a  variance 
estimate.  However  subseguent  experience  with  the  method 
indicates  that  it  probably  suffers  significantly  from  at 
least    two   factors: 

1.       the    planar   wavefront    assumption      probably    builds   in   a 
position    bias    as   in    the   case    of    the    L.S.    aethod;    and 
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2.  the  variance  estimate  is  inflated  since  part  of  the 
noise  being  measured  is  due  to  the  inadequacy  of  the 
planar   assumption. 

J.       COMPUTATIONAL    EXAHPLE 

for  a  quick  comparison,  the  three  methods  ieveloped  in 
this  chapter  are  now  applied  to  the  example  which  was  used 
in  Chapter  II.  .The  EXACT  results  calculated  previously  are 
included   in   Table  III  for   comparison. 


TABLE   III 

Single  Example  Comparison    of 
Plaiiar   Havefront   Methods 


TEa.D.§it   lijne  Elev.    Anqle 

Method  T   Tsecs)  S   TIea§f 

L.S.     •  0.67529319  15.22798 

L.S.C.  0.67527149  15.26372 

M.L.P.  0.67529007  15.10437 

EXACT  0.67527043  15.27002 

Apparent  Position        Estimate 

L.S.  (     1003.809  ,       3019.751     ,       866.1!  6    ) 

L.S.C.  (    1006.089  ,       3018-266    ,       868.235    ) 

M.L.P.  (       997.722  ,       3023.685    ,       859.355    ) 

EXACT  {    1005.966  ,       3017.899    ,       863.555    ) 


As  can  be  seen  easily,  the  M.L.P.  method  fares  rather 
poorly  in  all  regards,  even  worse  than  the  L.S.  method. 
This  will  be  confirmed  by  the  evaluations  made  in  Chapter  V. 
Also    worthy   of      note    is    the    apparent    tendency      o;     the    L.S.C. 
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method  to  correct  the  L.S.  method  back  toward  the  exact 
solution.  The  evaluations  of  Chapter  V  will  confirm  that  the 
L.S.C.  method  almost  always  yields  a  better  solution  than 
the    original  L.S.    solution. 
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IV.    A    SPHEEICAL    HAVEEIBQIT    MODEL 

A.       TEE    SPHEEICAL    WAYEFRONT    ASSUMPTION 

All  the  models  developed  in  Chapter  III  ire  limited 
primarily  by  the  assumption  that  the  wavefront  is  planar 
upon  arrival  at  the  hydrophone  array.  In  this  chapter  a 
model  is  developed  under  the  assumption  that  the  wavefront 
is      spherical   upon      arrival      at   the      array.  If    the      sound 

velocity  profile  were  constant  with  depth  then  the  spherical 
model  would  be  exact.  This  is  of  course  not  the  case,  but 
it  is  suspected  that  the  wavefront  is  better  modelled  as  a 
sphere  than  as  a  plane  because  that  small  piece  of  the  wave- 
front  which  passes  through  the  30  foot  cube  spanned  by  the 
array  is  locally  spherical.  That  is  because  et  ery  part  of 
that  piece  travelled  through  approximately  the  same  regions 
of    water,    experiencing   the    same   general   raybendiig    patterns. 

The  spherical  assumption  is  accurate  if  and  only  if  the 
speed  of  sound  is  constant  over  the  ray  path,  and  conseq- 
uently the  original  four  spherical  equations  apply  once 
again.      They  were: 


(4.1) 


(    X    -    D/2    )2    +    (    Y   +    D/2    }2    +    (    Z    +    D/2    )2      =      V2T2 

x 

(    X   +   D/2    )       +    (    y   -    D/2    )2    +    (    Z   +    D/2    )2      =      V2T2 

y 

(    X   +   D/2    )2    +    (    y   +    D/2    )2    +    (    Z    -    D/2    )2      =      V2T2 

z 

(    X   +    D/2    )2    +    (    y   +    D/2    )2    +    (    Z   +    D/2    )2      =      V2T2 

c 

It  has  been  stressed  previously  that  there  is  no  exact 
solution  [1,1, Z)  satisfying  all  four  equations  (4.1).  That 
is  because  the  time  values  on  the  right  hand  side  correspond 
to  the  reality  of  a  variable  velocity  profile.  However,  if 
the    spherical   wavefrcnt   assumption    is    to   be   accurate,    then  a 
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constant  velocity  is  the  assumed  case,  and  any  innacuracies 
in  the  time  values  are  regarded  as  due  to  timing  errors 
only.  Therefore,  under  the  spherical  assumption,  if  the 
true  time  values  Ui  (1=1,2,3,4)  were  known  and  substituted 
into  (4.1),  an  exact  solution  to  the  overdeter mined  system 
would  be  realized.  In  that  case  a.  solution  to  any  three  of 
the  equations  would  also  be  equal  to  that  unique  exact  solu- 
tion. In  particular,  the  NAVY  solution  of  Chapter  II  would 
be  the  true  solution.  Therefore  in  terms  of  the  coordinate 
system  referenced  to  the  c  hydrophone,  X  would  be  given  by 
(4.2). 

D  V2  ,       2  2 


+ 


(    U4    -   U2    )  (4.2) 


2  2D 

However,  X  is  also  given  by  (4.3),  whers  C1  is  the 
direction  cosine  along  the  X  axis  of  the  vector  from  the  c 
hydrophone    to    the   sound    source. 

x  =  v  u4  ci  •  («-3) 

The  time  value  of  interest  at  the  moment  is  J  4,  the  time 
at  the  c  hydrophone.  Therefore  let  U  =  U4  for  clarity  of 
presentation,  and  then  eguate  the  expressions  in  (4.2)  and 
(4.3)  in  order  to  solve  for  U.  If  this  same  lDgic  is  also 
applied  to  the  similar  expressions  for  the  distances  to  the 
y   and    z    hydrophones,    the   results    are    (4.4). 


U.       =     y  IT    -     (2DUC./V)    +     (D/V)2  i   =    1,2,3  (4.4) 

The    expression     (4.4)       will   be      useful    in    the    development 
of    the    model  of    this   chapter. 


B.       LEAST    SQUARES    MODELS 

A    direct      approach   might    be      to    apply    the      Is ast   squared 
error    technique      to    the      spherical   equations,         in    a      manner 
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paralleling  that  used  on  the  four  planar  equations  of  the 
L.S.  model  in  Chapter  III.  However  the  formulae  and  equa- 
tions that  result  are  exceedingly  complex,  involve  fourth 
degree  powers  of  the  data  values  Ti,  and  have  thus  far 
defied  all  solution  attempts.  Therefore  this  ids  a  was  aban- 
doned in  favor  of  the  maximum  likelihood  approach  which 
follows. 

C.   MAXIHOB  LIKELIHOOD  COMPUTATIONS 

As  in  the  M.L.P.   model,   Gaussian  noise  is  assumed  for 
the  time  data  values.   Therefore 

Ti   =   Ui   +   Ei        i  =  1,2,3 
where  the  Ei  are   independent  identically  distributed  normal 
random   variables  with  mean   zero   and  variance   S2.     The 
density  of  each  Ti  is  therefore  (4.5)  . 


f±(T)   = 


—  exD 


27T  s 


~    (  T   -  U   )2 
2  S      X     X 


i  =  1,2,3,4    (4.5) 


To      form  the      likelihood      function,         multiply    the      four 
densities    together,    to  obtain     (4.6). 

Form   the  log   likelihood    function   by   taking    nitural    loga- 
rithms   of     (4.6)     to    obtain   L1    of     (4.7)  . 

r    3 


L        =      -2    ln(27T)    -    4    ln(S)    - 


2    S* 


Z2  2 

(T.-U.)       +     (1\-U) 
11  4 

.=  1 


(4.7) 


In  order  to  maximize  L1  with  respect  to  C1,  C2,  C3  and 
U,  it  is  sufficient  to  minimize  L2  of  (4.8),  where  a  substi- 
tution  for    Ui   has   been  made    using    (4.4). 
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3     r 


L_      =       (T   -U) 
2  4 


Y]    T.    -    yv"   -    (2DUC. 


1=1 


/V)     +     (D/V) 


(4.8) 


Now    add    the    usual   direction      cosines    constraint    and    form 
the    lagrangian    function   I   of    (4.9). 

p  3   _ 


-x 


iti x 


(4.9) 


For  notational  convenience,  let  Ki  be  given  by  (4.10) 


K 


i   =   U   -  (2  D  U  C\  /  V)  +  (D/V)2       i  =  1,2,3     (4.10) 

Take  the  partial      derivative    of    L    with    respect      to    Ci    to 
get     (4.11). 


3l 


f2    (T.  -   K.) 

l      l 


-2   K. 

l 


-2  D  U 


-  2  \  C.  (4.  11) 

A  i     i  =  1,2,3    V     ' 


Simplify  (4.11)  and  equate  to  zero  to  yield  (4.12). 

d  u  a  -y*i 


c.  = 


xv  vv*r 


i  =  1,2,3 


(4.12) 


Multiply  the  three  equation  in  (4.12)   by   Ci  (1=1,2,3) 
respectively,    and   add  them    together.     Than   use   the 
constraint  on  the  sum  of  squares  of  the  direction  cosines  to 
obtain  (4.13). 

Substitute  (4.13)  into  (4.12)  to  get  the  maxi  mum  likeli- 
hood estimate  for  the  direction  cosines,  as  in  (4.14). 
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Ti    -VKx 


(4. 14) 


V^X 


T:  'VS 


i    =    1,2,3 


Now    take  the    partial  derivative    of      L    with    r? spect    to    U, 
as    in    (4.15). 

3 


3l 


i=l  u 


r/Td  -V^ 


(2U    -    2DC./V) 


+    2     (U    -   T    )       (4«15) 


Equate  (4.15)  to  zero  and  solve  for  0  to  obtain  the 
maximum  likelihood  estimate  of  the  time  value  (4.16)  at  the 
c  hydrophone. 

3 


u  = 


•  *  M* 


(4.16) 


T.  -n/k" 
1       V  D 


*  -     I 

j  =  l 

Finally  take  the  partial  derivative  of  11  in  (4.7)  with 
respect  to  S.  Equate  it  to  zero  and  solve  to  get  (4.17), 
the  maximum  likelihood  estimate  of  the  variance. 


r  3 


Y   (T.  -   K.  )    +   (T   -  U) 


Lj=i 


(4.17) 


D.       SOLUTION   BY    A    MODIFICATION    OF    NEWTON'S    METHOD 

The  solutions  given  by  equations  (4.14),  (4.16),  and 
(4.17)  once  again  form  a  set  of  equations  which  would  seem 
to  be  solvable  by  natural  iteration  as  in  the  M.^.P.  method. 
Unfortunately  this  time  the  technique  fails  to  converge.  A 
mathematical  tool  is  needed  which  is  stronger  t han  natural 
iteration.         What      is   used      is   a      modified    four      dimensional 
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version   of    Newton's      Method    [Eef.    5    p. 47]    to      sei rch    for    the 
roots   of   a   set    of   four  equations. 

The    objective   is    to    determine   the    values    C1,    C2,      C3    and 
0   which    satisfy     (4.14)      and     (4.16).         For   further    notational 
convenience  define      the   values      M   and      N   as      in    (4.18)         and 
(4.19). 


M 


(4.  18) 


(4.19) 


Given   those    definitions    of    M      and    N,       then    tie    equations 
whose    roots  are    desired   can    be   simplified    to    (4.20). 


u 


=    h.(crc2,c3,u) 


=    h4(crc2,c3,u) 


i  -V*T 


v^r 


T4    - 


i   =    1,2,3 


(4.20) 


D   M 
V 


1    -    N 


Now  define  error  functions  as  in  (4.21).  Tha se  evaluate 
the  amount  of  error  in  each  of  the  equations  (4.20)  for  any 
set    of   values   for    (C  1,C2,C3 , U) . 


gi 


=     c.   -  h. 
i        i 


i   =    1,2,3 


Let        G        be 

(g1,g2,g3,g4)  <  . 


(4.21) 

^4         =      U    "    h4 

the     four        dimensional        coli mn        vector 
Finally      let    GP    be      the    matrix      of    partial 
derivatives   of    G,    as   in    (4.22). 

Newton's  Method  in  four  dimensions  says  that  if  X(n)  J-s 
a  four  dimensional  vector  holding  the  current  approximate 
roots  C1,C2,C3  and  Ur  then  X(n+1)  will  be  an  improved 
answer,    where    X(n+1)     is    given    by     (4-23). 


53 


GP    (C    ,C 2'c3'u> 


3gl 

a*3 

3gl 

Bgi 

3C1 

^S 

^C3 

3C4 

9^2 

dc2 

Bc2 

eg2 

Bc3 

Sg2 
^C4 

Bg3 

3g3 

Bg3 

Bg3 

3C1 

Sc2 

3C3 

SC4 

Bg4 

3g4 

Bg4 

3g4 

dc 


d=. 


(4.  22) 


-n+1 


X 
— n 


-    [GP<^»] 


-1 


(V 


(4.23) 


Of  course  it  is  necessary  to  calculate  the  derivatives 
held  in  the  matrix  GP  in  order  to  use  (4.23).  Ti  ose  deriva- 
tives  are    given    in    Appendix    B. 

Unfortunately  when  multidimensional  versions  of  Newton's 
Method  are  applied,  there  is  often  a  tendency  for  the  method 
to  converge  slowly,  or  even  diverge.  This  is  because  it 
tends  to  overshoot  the  best  answer  for  each  iteration.  To 
alleviate  this  problem,  a  modification  is  made  to  the 
method.  At  the  end  of  each  Newton  iteration,  prior  to 
proceeding  with  the  next  iteration,  a  Golden  Section  Search 
[Ref.  6  ]  is  performed  to  find  the  best  possible  answer  in 
the  direction  of  the  new  iterative  solution.  Specifically 
the  line  in  four  dimensional  space  from  X(n-1)  of  the 
previous  iteration  to  X(n)  of  the  present  i:  eration  is 
searched  for  the  best  answer.  The  definition  o£  the  'best* 
answer  is  that  point  along  the  search  line  whi^ h  minimizes 
the    sum   cf    squared    error   functions,    namely    (4.24| . 


f   2 


(4.  24) 
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The  current  iterative  solution  X(n)  is  thea  given  the 
value  of  the  minimizing  point  resulting  from  the  Gclden 
Section  Search.  Then  the  next  iteration  of  New:on's  Method 
is  performed,  along  with  another  Golden  Section  Search  to 
find    the  next  iterative    solution  X(n+1). 

E.       THE    MAXIMUM    1IKE1IHOOD    METHOD 

In  summary,  the  fourth  and  final  alternative  method  to 
estimate  apparent  positions    is   as   follows: 

1.  let    0  =    T4    initially; 

2.  use  the  L.S.  method  (Se321)  to  initially  estimate  the 
values  of   Cj     (j=1,2,3); 

3.  set    X(1)       =        (01,02,03,0)  ; 

4.  initialize  the  Newton  Iteration  counter  :  I  =  1  ; 

5.  calculate  the  values  Ki,  M  and  N  in  accordance  with 
equations  (4.10),  (4.18)  and  (4.19); 

6.  calculate  the  error  function  vector   G  (X( I) )  ,   using 

(4.20)  and  (4.21)  ; 

7.  calculate  the  derivative  matrix  GP(X(I)),  using  the 
results    in  Appendix    B; 

8.  invert   GP  (X  (I))  ; 

9.  calculate  the    new   Newton   estimate   X(I+1)     from    (4.23); 

10.  perform  a  Gclden  Section  Search  along  the  line 
between  X(I)  and  X(I  +  1)  to  find  the  point  which  mini- 
mizes   (4. 24)  ; 

11.  let  X(I  +  1)  be  equal  to  the  minimizing  point  found  in 
the    previous    step; 

12.  increase    the    Newton    iteration    counter    :    I    =   1+1     ; 

13.  reiterate  steps   5   through    12    until    the    values 

X(I)       =       (01, 02, 03,0) 
converge    within   acceptable   tolerances; 

14.  calcualte    the    estimate   of    S2    using    (4.17).. 
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This  method  shall  hereafter  be  referred  to  as  the 
'Maximum  Likelihood  Spherical  Method',  or  M.L.S.  for  short. 
As  will  te  seen  from  the  evaluations  made  in  Chapter  V,  the 
M.L.S.  method  is  apparently  the  only  alternative  to  consis- 
tently rival  the  performance  of  the  currently  used  NAVY_A 
method.  It  has  the  additional  advantage  that  it  estimates 
the    error   present  in    the   time   data   values. 

F.       A    COHPOTATIONAL    EXAMPLE 

This  latest  method,  M.L.S.,  is  now  applied  to  the  same 
example  considered  in  Chapters  II  and  III.  'or  a  quick 
comparison,  Table  IV  lists  the  results  of  using  all  the 
methods.  The  error  vector  lengths  are  the  distances  of  each 
position  estimate   from   the    EXACT    apparent    position. 

Since  this  is  cnly  one  example,  this  title  is  not 
presented  for  the  purpose  of  any  broad  conclusions.  However 
it   is    of   interest   to   note   that    in    this   example 

1.  the    M.L.S.      method   outperforms   all    others,      including 
the    NAVY_A  method;    and 

2.  in   many    ways    the   original   NAVY    method   outperforms   the 
adjusted    NAVY_A  method. 

G-       VARIANCE   ESTIMATICH 

In  the  computational  example,  the  time  data  values  used 
were  exact  since  the  methods  of  Appendix  A  coi Id  be  used 
with  the  known  linear  velocity  profile.  Therefore  the 
appropriate  value  for  variance  in  the  time  vali es  would  be 
zero.  For  this  error  free  example,  the  estimates  shown  in 
Table  V  were  obtained  for  the  standard  deviations  of  timing 
noise,    using   the    two    maximum   likelihood    methods. 

In  the  case  of  this  one  example,  the  spherL  cal  assump- 
tion is  apparently  an  improvement  over  the  planar,  since  the 
M.L.S.    estimate    of    error   is    only    4%   of   the    M.I. P.  .     estimate. 
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TABLE    IV 

Single   Example   Comparison   of   All 

Methofis 

Method 

Transit  Time 
T  "decs) 

Elev.    Anc[J 
1    C2§2sf 

.e 

Length    of 
l££2f    Vector 

NAVY 

0. 67526969 

15.26459 

0.412  8 

NAVY_A 

0. 67527027 

15.26461 

0.4840 

L.  S. 

0.67529319 

15.22798 

3.739 

L. S.C. 

0.67527149 

15.26372 

0.522 

a. i.  p. 

0. 67529007 

15. 10437 

13.64) 

M.L.S. 

0. 67527C49 

15.26677 

0.411 7 

EXACT 

0. 67527043 

15.27002 

A  pparent 

Position 

Estimate 

NAVY 

(    1005.957 

,      3017.8  74 

/ 

868.14  3    ) 

NAVY_A 

(    1006.087 

,       3018.259 

/ 

868.25  5    ) 

L.S. 

(    1003.809 

,       3019.751 

/ 

866.126    ) 

L  •  S.  C  ■ 

(     1006.089 

,       3018.266 

/ 

868.235    ) 

M.I.?. 

(       997.722 

,       3023-685 

/ 

859.35  5    ) 

M.L.S. 

(    1006.195 

,      3018.190 

$ 

868.375    ) 

EXACT 

(    1005.966 

,       3017.899 

i 

868.55  5    ) 

That  is  regarded  as  an  improvement  because  the  correct 
answer  is  zero.  The  higher  M. L- P.  estimate  is  indicative  of 
the    inflation   due    to    the   planar   assumption. 

The  tracking  range  which  provided  data  for  this  study 
records  time  values  to  seven  decimal  places.  Tierefore  the 
standard  deviation  estimated  by  the  M.L.S.  method  is  of 
particular  interest,  since  it  indicates  errors  in  the 
seventh  decimal  place  even  when  there  is  no  error  present. 
Since    the   data    was    actually    error    free    in    this   es ample,       the 
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TABLE    V 

Maximum    likelihood   Error    Estimates 
for  the  Example   Problem 


Model 

Planar 

Spherical 


Method 

M.L.P. 
M.L-S. 


Est.    Std.    Devil  tion 
5-517    E-6        sees. 
2.191    E-7        se:s. 


estimate  is  a  measure  of  the  variation  induced  b{  the  spher- 
ical vavefront  assumption.  A  broader  discussion  of  error 
estimation    will    be    undertaken   in    Chapter    VI. 
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V.    EVALUATIONS    OF    MODELS 

A.  GENEEAL 

There  are  many  sources  of  errors  in  the  overall  hydro- 
phonic   tracking    problem.      These   include,    among   others: 

1.  errors  in   estimation    of   the   sound   velocity    profile; 

2.  inhomogeneity  of  the  velocity  profile  over  time  and 
horizontal   displacement;    and 

3.  possible  errors  in  measuring  the  positiois,  and  the 
angles  of  tilt  and  rotation  for  the  hydrophonic 
arrays. 

This  study  focuses  on  those  errors  which  :>ccur  during 
the  computations  preceding  the  ray  tracing  procedure.  To 
evaluate  the  performance  of  the  methods  developed  in 
Chapters  II,  III  and  IV,  it  is  necessary  to  control  strictly 
the  ray  tracing  procedure.  Only  in  that  way  can  the  differ- 
ences found  between  methods  be  attributed  to  the  differences 
between  models,  and  not  to  any  source  of  error  oi tside  those 
methods. 

It  was  originally  hoped  that  the  various  methods  might 
be  compared  by  applying  them  to  real  tracking  data.  However 
it  was  found  that  the  overall  tracking  problem  a  ad  too  many 
large  sources  of  error  to  allow  the  methods  to  demonstrate 
any  differences.  Therefore  the  methods  were  coa  pared  under 
a   more    tightly    controlled  simulated   environment. 

B.  SIBUIATION    SCENABIO 

Two  different  simulations  are  used  to  compare  the  six 
different    methods.  Both    use    a      basic    scenario      similar   to 

that    of    the      computational    example    explored    in      Chapters   II, 
III    and   IV.      That   example   assumed    that: 
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1.  the  velocity  versus  depth  relationship  is  linear,  and 
given  exactly   by: 

V       =      4840.7      +       0.03314    •    DEPT1     ; 

2.  the  acoustic  center  of  each  array  is  at  a  depth  of 
1300    feet; 

3.  the  hydrophone  arrays  are  ail  level,  and  their  X,  Y 
and  Z  arms  are  parallel  to  the  respective  coordinate 
axes   of    the    tracking    range. 

Under  these  circumstances  the  methods  se:  forth  in 
Appendix  A  can  be  used  to  compute  the  exact  values  for  the 
hydrophone  times,  ray  transit  time  and  elevation  angle  for  a 
sound  wave  emanating  from  a  source  at  any  specified  loca- 
tion. Therefore  when  the  methods  are  applied  to  those  exact 
times,  the  resulting  estimated  positions  can  be  compared  to 
the    known   true    position. 

C-       SIMULATED    EBROR    FOE    TIMING    DATA 

The  models  developed  in  this  study  were  designed  to 
improve  position  estimation,  especially  in  the  presence  of 
errors  in  the  timing  data.  Lacking  any  better  hid  del  at  this 
time  for  timing  errors,  the  simulated  environment  includes 
an  assumption  of  Gaussian  errors  for  the  hydrophone  times. 
Therefore  realistic  timing  data  can  be  simulate!  by  adding 
to  each  exact  time  value  a  random  quantity  of  normally 
distributed  error.  The  mean  of  the  error  is  assumed  to  be 
zero.  The  variance  was  estimated  from  real  tracking  data, 
using  the  variance  estimating  property  of  the  H.L.S.  model. 
The  data  from  one  tracking  run  was  used,  involving  six 
hydrophone  arrays  and  733  position  estimates  (sea  Chapter  VI 
for      data      selection      details) .  Each      position      from      the 

tracking  run  produces  one  estimate  of  the  variance.  The 
variance  value  chosen  for  use  in  the  simulations  was  the 
median  of  the  733  variance  estimates  produced  bf  the  M.L.S. 
method.       That   value    was 
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S2         =    9. 1204   E-12    secs2  . 
That  is  the  same  as  a  standarad  deviation  of 

S     =   3. 02   E-6    sees  . 
Each  of  the  two  types  of   simulations  was  rui  four  sepa- 
rate times.    Each   run  was  done  with   a  different  specified 
variation.    Those  four  distinct  error  conditions  are  delin- 
eated in  Table  VI  . 


r 

TABLE    71 
Simulation   Error 

Levels 

EUN 

LEVEL 

VARIANCE 

STD. 

DEVIAC ION 

1 

zero 

0.0 

3.0 

2 

low 

9.12    E-14 

3.02 

E-7 

3 

medium 

9.  12    E-12 

3.02 

E-6 

a 

high 

9.1 2    E-10 

3.02 

E-5 

D.   SINGLE  AERAY  SIMULATION 

In  the  first  simulation,  the  intent  was  to  compare  the 
methods  pairwise,  so  as  to  determine  which  me:  hod  is  more 
likely  to  produce  the  more  accurate  estimate  of  a  given 
sound  source  position.  One  thousand  positions  were  chosen 
at  random.  Each  position  was  3000  feet  from  the  array.  The 
positions  were  uniformly  spread  over  the  surface  of  a  sphere 
of  radius  3000  feet,  centered  at  the  array,  truncated  above 
by  the  water  surface  (depth  0)  and  below  by  the  depth  of  the 
array  (1300  feet).  The  methods  set  forth  in  Appendix  C  were 
used   to  assure   that   the   random  positions   selected   were 
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uniformly  distributed  over  the  surface  area  of  the  truncated 
hemisphere. 

Each  of  the  1000  randomly  chosen  source  pDsitions  was 
then  processed  as  follows: 

1.  calculate  the  exact  hydrophone  times,  using  the 
methods  of  Appendix  A; 

2.  add  to  each  of  the  four  exact  times  a  random  value  of 
error  at  the  specified  level  (zero,  low,  medium  or 
high) ; 

3.  apply  each  of  the  six  methods  to  the  hydrophone 
times,  generating  six  different  apparent  positions; 

4.  apply  the  ray  tracing  procedure  to  each  of  the  six 
positions,  using  layers  that  are  25  feet  thick,  and 
utilizing  the  known  linear  velocity  profile,  thereby 
producing  six  different  estimates  of  the  sound  source 
location ; 

5.  compare  the  six  different  estimates  pairwise  to  see 
which  method  in  each  pair  produced  the  estimate 
closest  to  the  true  sound  source  location. 

The  comparison  being  made  is  that  one  method  is  consid- 
ered preferrable  to  the  other  if  it  more  frequently  produces 
the  more  accurate  estimate. 

The  layer  thickness  of  25  feet  was  selected  order  to 
simulate  the  actual  procedure  at  the  tracking  range  which 
provided  data  for  this  study.  However,  as  discussed  in 
Chapter  I,  when  the  velocity  profile  is  smooth  and  known 
exactly,  the  process  is  very  robust  with  respect  to  the 
thickness  used.  The  thickness  values  1,  10,  and  25  were  all 
tried,  with  virtually  no  changes  in  any  of  the  comparisons 
between  methods. 
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E.       SINGLE    ARRAY    SIMULATION    RESULTS 

Tables  VII  and  VIII  contain  the  results  of  the  single 
array  simulation.  Each  tabled  entry  represents  :  he  fraction 
of  time  that  the  method  of  that  row  produced  a  better  esti- 
mate than  the  method  of  that  column.  For  example,  in  Table 
VII,  the  l.S.  method  outperformed  the  M.L.S.  method  in  only 
5.8%    of    the    1000    trials    with   low    error    values. 

In  a  one  tailed  test  that  one  method  is  setter  than 
another,  these  binomial  proportions  are  significant  at  the 
0.0  5  level  if  they  exceed  0.526.  Symmetrically,  one  method 
is  significantly  worse  than  another  at  the  0.05  level  if  the 
proportion      is   less      than   0.474.  For    the      0.01     level      the 

corresponding  critical  values  are  0.537  and  0.463 
respectively. 

The   results    indicate   that: 

1.  as  previously  claimed,  the  NAVY_A  method  usually 
outperforms  the  unadjusted  NAVY  method;  of  particular 
interest  is  the  case  of  zero  error  whi; h  actually 
compares  the  relative  ability  of  each  method  to 
produce  the  exact  answer  when  given  the  exact  times; 
in  those  cases  the  NAVY_A  method  does  extremely  well 
against    the    NAVY    method; 

2.  under  all  error  conditions  the  most  successful 
performer  is  the  M.L.S.  method,  since  it  always  has  a 
favorable  (greater  than  0.5)  comparisDn  fraction 
against  all  other  methods;  the  M.L.S.  fractions  vary 
little   over    the  four    error    levels; 

3.  under  all  error  conditions,  the  worst  performer  is 
always   the    M.L.P.    method; 

4.  spherical  methods  consistently  outperform  planar 
methods;    and 

5.  increased  error  levels  tend  to  lessen  the  distinction 
between    methods;      in    the      zero   error   case    comparisons 
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TABLE    VII 

Single    Array    Simulation    Results 
for  Lower   Error  Levels 

ERROR 

LEVEL    : 

ZERO 

NAVY 

NAVY_A      I 

,.S.            I 

j  •  O  «   \*  • 

M.L.P. 

NAVY 

.  011 

.950 

.676 

.987 

NAVY_A 

.98  9 

.950 

.728 

.987 

L.S. 

.05  0 

.050 

.050 

.987 

L.S.C. 

.32  4 

.  272 

.950 

.987 

M.I. P. 

.013 

.  013 

.013 

.013 

M.L.S. 

.580 

.566 

.943 

.654 

.987 

EEROR 

LEVEL    : 

iow 

NAVY 

NAVY_A      I 

,.S.            I 

j  •  O  •  C  « 

M.L.P. 

NAVY 

.  165 

.952 

.648 

.987 

NAVY_A 

.835 

.952 

.693 

.987 

L.S. 

.048 

.048 

.048 

.987 

J-t    *    *-  •    v-   * 

.352 

.307 

.952 

.987 

M.I. P. 

.013 

.  013 

.013 

.013 

M.I.S. 

.604 

.595 

.942 

.631 

.987 

M.L.S. 
.420 
.434 
.057 
.346 
.013 


M.L.S, 
.396 
.405 
.058 
.369 
.013 


are  in  the  interval  (0.01,0.99),  while  in  the  high 
error  case  that  interval  is  narrowed  considerably  to 
(0.37,0.63)  . 

F.       DOUBLE    ABRAY    SIflOIATION 

In  the  second  simulation,  the  intent  again  was  to 
compare  the  methods  pairwise,  this  time  determining  which 
method  is  mere  likely  to  produce  positions  which  agree  more 
closely   in    the    two    array   crossover   problem.      This   is   not   the 
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1 

TABLE    VIII 

Si 

.ngle    Array    Simulation    Results 
for  Higher  Error   Levels 

ERROR 

LEVEL 

:       MEDIUM 

NAVY 

NAVY_A      L.S. 

L.S.C. 

M.L.P. 

NAVY 

.464 

.819 

.559 

.987 

NAVY_A 

.546 

.821 

.566 

.987 

L.  S. 

.181 

.  179 

.181 

.971 

L.  S.C. 

.441 

.  434 

.819 

.987 

H.I. P. 

.013 

.013 

.029 

.013 

H.I.S. 

.56  5 

.  563 

.823 

.550 

.987 

EEEOR 

LEVEL 

:       HIGH 

NAVY 

NAVY_A      L.S. 

L.S.C. 

M.L.P. 

NAVY 

.  486 

.544 

.439 

.562 

NAVY_A 

.514 

.546 

.516 

.559 

±l  a    S. 

.45  6 

.454 

.457 

.557 

L.S.C. 

.56  1 

.484 

.543 

.5  62 

M.I. P. 

.438 

.  441 

."44  3 

.438 

H.I.S. 

.56  9 

.569 

.600 

.573 

.627 

M.L.S 
.435 
.  437 
.  177 
.450 
.013 


M.L.S 
.431 
.431 
.400 
.427 
.373 


same  question  as  that  addressed  by  the  single  array  simula- 
tion. The  two  estimates  produced  by  any  one  method  may  be 
very  close  to  each  other,  and  yet  be  far  away  from  the  true 
position. 

For  the  double  array  scenario  two  arrays  are  used,  sepa- 
rated by  7500  feet,  toth  at  depths  of  1300  feet..  The  arms 
of  each  array  are  parallel  to  the  corresponding  coordinate 
system  axes.  Once  again  10  00  positions  were  ranlomly  gener- 
ated   in    a   uniform   manner,    this    time    over    a    3    dimensional   box 
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Figure    5.1        Double   Array   Simulation   Configuration. 

running  crossways  between  the  two  arrays.  The  box  is  1300 
feet  deep,  5000  feet  long  and  1000  feet  across  (see  figure 
5.1)  . 

Each   of   the    1000    randomly      chosen    sound    source    positions 
in    the    two    array   simulation    were    processed    as   foLlows: 

1.      calculate   the      exact    hydrophone      times   for      the    first 
array  using    the  methods   set   forth   in   Appendix    A; 


66 


2.  add  to  the  exact  times  some  random  error  at  the  spec- 
ified level; 

3.  Repeat  steps  1  and  2  for  the  second  array ,  usin^  a 
new    set    of   random   error   values; 

4.  apply  each  of  the  six  methods  to  the  time  values  from 
both  arrays,  producing  six  differen:  pairs  of 
apparent    positions; 

5.  apply  the  ray  tracing  procedure  to  bo: h  apparent 
positions  in  each  of  the  six  pairs,  utilizing  the 
known  linear  velocity  profile,  producing  six  pairs  of 
estimated   sound   source    positions; 

6.  for  each  of  the  six  pairs  of  positions,  cilculate  the 
distance    between    the    two    positions    in    the    pair; 

7.  make  pairwise  comparisons  of  the  six  different 
distances,  to  see  which  method  in  each  pair  exhibits 
the  closest  agreement  between  its  two  position  esti- 
mates. 

This  time  the  comparison  being  made  is  that  one  method 
is  considered  -better  than  another  if  the  positions  it 
produces  agree  more  closely  more  often  than  those  of  the 
other    method. 

G.       DOUBLE    ARRAY    SIMD1ATI0N    RESULTS 

Tables  IX  and  X  contain  the  results  of  the  1 ouble  array 
simulation.  Each  tabled  entry  represents  the  fraction  of 
time  that  the  method  of  that  row  produced  a  pair  of  esti- 
mates which  were  in  closer  agreement  than  the  estimates 
produced      by      the      method      of      that      column.  Significance 

criteria  for  these  proportions  are  the  same  as  for  the 
single   array  simulation. 

These  results  are  similar  to  those  of  the  single  array 
simulation,    because    they   indicate    that: 
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TABLE   . 

IX 

Double    Array    Simulat 
for  Lower  Error 

ion    Results 
Levels 

ERE  OP. 

LEVEL 

I 

ZERO 

NAVY 

NAVY_A 

L.S. 

L.S.C. 

M.L.P. 

M.L.S. 

NAVY 

.  366 

.989 

.726 

.968 

.212 

NAVY_A 

.634 

.989 

.734 

.968 

.212 

L.  S. 

.01  1 

.  011 

.011 

.555 

.014 

L.S.C. 

.27  4 

.  266 

.989 

.968 

.  193 

M-L.P. 

.03  2 

.  032 

.445 

.032 

.033 

M.L.S. 

.78  8 

.788 

.986 

.807 

.967 

ERROR 

LEVEL 

* 

LOW 

NAVY 

NAVY_A 

L.S. 

L.S.C. 

M.L.P. 

M.L.S. 

NAVY 

.448 

.989 

.661 

.952 

.310 

NAVY_A 

.55  2 

.989 

.684 

.952 

.  312 

L.S.. 

.01  1 

.  011 

.011 

.560 

.014 

L.S.C. 

.33  9 

.316 

.989 

.952 

.254 

M.I. P. 

.04  8 

.048 

.440 

.048 

.038 

M.I.S. 

.690 

.  688 

.986 

.746 

.962 

1.  the  NAVY^A  method  outperforms  the  original  unadjusted 
NAVY  method,  although  the  difference  is  not  signifi- 
cant  in    the    higher   error   level    cases; 

2.  the  most  successful  performer  is  consistently  the 
M.L.S.    method; 

3.  spherical  methods  almost  always  outperform  planar 
methods; 

4.  increased  error  levels  tend  to  lessen  the  distinction 
between    methods. 
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TA3LE    X 

Double    Array    Simulation    Results 
for  Higher  Error   Levels 


ERE  OR 

LEVEL 

:       MEDIUM 

NAVY 

NAVY_A 

L.S. 

L.S.C. 

M.L.P. 

M.L.S 

NAVY 

.498 

.841 

.513 

.808 

.443 

NAVY_A 

.502 

.842 

.542 

.809 

.440 

L.S. 

.159 

.  158 

.159 

.525 

.  160 

L.S.C. 

.437 

.458 

.84  1 

.808 

.433 

M.I. P. 

.192 

.  191 

.475 

.192 

.  144 

H-I.S. 

.557 

.560 

.840 

.567 

.856 

NAVY 

NAVY_A 

.513 

L.S. 

.444 

L.S.C. 

.465 

M.L.P. 

.521 

M.L.S. 

.55  8 

ERROR    LEVEL  ;       HIGH 

NAVY          NAVY_A      L.S.  L.S.C.  M.L.P.  M.L.S 

.487           .556  .535  .479  .442 

.558  .499  .481  .440 

.442  .445  .451  .422 

.501            .555  .480  .438 

.519           .549  .520  .428 

.560            .578  .562  .572 


There  are  however  some  indications  from  these  results 
which  are  different  from  those  of  the  single  array  simula- 
tion,   such    as: 

1.  the  worst  performer  was  consistently  the  L.S.  method, 
rather   than    the    M-L.  P.    method; 

2.  in  the  high  error  case  the  M.L.P.  method  is  equal  to, 
or  even  marginally  better  than,  every  other  method 
except    the    M.L.S.    method; 
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3.  the  performance  of  the  M.L.S.  method  is  noticeably 
better  under  error  free  conditions. 
The  last  two  inferences  are  perhaps  the  most  inter- 
esting. The  maximum  likelihood  approach  was  intended  to 
estimate  and  account  for  Gaussian  errors  in  the  timing  data 
values.  Hence  it  is  really  not  very  surprisii g  that  the 
M-L.S.  method  does  well  with  error  prone  data.  But  it  is 
interesting  to  note  that  the  M.L.P.  method  also  does  well 
under  high  error  levels,  even  though  it  probably  suffers 
from  a  bias  due  to  the  planar  assumption. 

Of  even  greater  interest  is  that  the  M.L.S.  method  seems 
to  be  at  its  best  when  compared  to  other  methods  under  error 
free  conditions.  This  result  was  unexpected,  and  indicates 
that  the  M.L.S.  method  not  only  handles  timing  errors  well, 
as  was  intended,  but  apparently  also  does  an  evei  better  job 
of  approximating  the  elusive  exact  solution  to  the  original 
four  spherical  eguations  of  (4.1)  and  (2.1)  when  the  exact 
time  values  are  available. 

H.   LIMITATIONS  ON  INTERPRETATION  OF  RESULTS 

The  results  of  both  simulations  seem  to  imply  that  the 
M.L.S.  method  outperforms  the  currently  used  NAVY_A  method 
on  any  randomly  chosen  sound  source  position,  with  or 
without  timing  errors.  These  are  encouraging  results. 
Nevertheless  the  reader  is  cautioned  that  thess  tests  are 
just  simulations,  and  like  all  simulations,  must  make 
assumptions  which  cannot  fully  reflect  the  reali:  y  of  actual 
hydrophonic  tracking  conditions.  The  most  important  assump- 
tions made  for  these  simulations  are: 

1.  the  sound  velocity  profile  is  known  exactly,  and  is 
linear; 

2.  the  errors  in  the  timing  values  from  any  hydrophone 
are  normally  distributed  with  mean  zero,  and  are 
independent  of  the  noise  in  any  other  hydrophone;  and 
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3.   the   parameters  of   the  simulation   were  fixed;    for 

example   the  single   array   simulation   used  a   fixed 

range  of  3000  feet,   and  both  simulations  used  arrays 

at  1300  feet  of  depth,  with  fixed  orientations  to  the 

range  coordinate  system. 

The  first  assumption  is   probably  of  little  consequence. 

It  not  only  greatly  facilitates  computations,  but  also  helps 

to  isolate   the  initial  angle  and   time  estimati. on   problem 

from  the  unrelated   errors  involved  in  the   velocity  profile 

estimation  procedure. 

The  second  assumption  is  somewhat  more  troublesome. 
Errors  may  not  be  Gaussian  at  all,  or  if  they  are,  the  mean 
may  not  be  zero.  Unfortunately  each  positioi  estimation 
involved  only  four  equations,  and  therefore  did  r.ot  allow 
for  estimation  of  more  than  four  paramenters.  Tierefore  the 
error  mean,  being  a  fifth  parameter,  could  not  be  estimated. 
Also  the  errors  of  any  one  hydrophone  may  very  well  not  be 
independent  of  the  errors  of  the  other  three  phones  on  that 
array.  Fortunately  these  concerns  are  offset  somewhat  by 
the  results  of  the  double  array  simulation,  whs  rein  it  was 
found  that  the  M.L.S.  method  was  at  its  best  when  there  was 
no  noise  at  all. 

Also  the  error  type  and  level  may  depenl  on  other 
factors,  such  as  the  target's  range,  elevation  and  azimuth 
angles  from  the  array.  This  highlights  the  concerns  of  the 
third  assumption.  There  is  considerable  room  here  for 
future  work  concerning  the  dependency  of  resuLts  on  such 
complicating  factors. 

Lastly  it  should  be  pointed  out  that  the  simulations 
make  comparisons  only  on  the  binomial  basis  of  better  versus 
worse  in  1000  trials.  The  magnitudes  of  the  actual  differ- 
ences are  ignored.  It  is  possible,  though  perhaps  unlikely, 
that  while  one  method  marginally  outperforms  a  second  method 
in  most  trials,  in  all  the  remaining  trials  the  :  irst  method 
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is  much  worse  than  the  second 
further  work. 


There  is  also  room  here  for 
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VI.     CONCLUSIONS    AND    RECOMMENDATIONS 

A.       ESTIMATION    OF   TIMING    DATA    VARIATION 

One  of  the  primary  purposes  of  this  study  #as  to  esti- 
mate the  amount  of  variability  in  the  timing  data  teing 
recorded      during   actual      tracking   runs.  This    problem      was 

addressed  by  the  M.L. P.  and  M.I.S.  maximum  likelihood 
models.  The  M.L. P.  model,  as  previously  discussed,  suffered 
from  a  bias  due  to  the  planar  assumption,  was  the  poorest 
estimator  of  positions  among  all  the  models,  and  produced  an 
inflated  variance  estimate.  Therefore  the  spherical  model 
M.L.S.      is    used    to    estimate    the   data   variance. 

The  variability  that  is  measured  by  the  M.L.S.  model  is 
made  up  of  three  components.  First  there  is  :  he  variance 
induced  by  the  spherical  assumption.  Then  there  is  the 
variance  caused  by  the  seven  decimal  accuracy  used  when 
recording  the  data.  Finally  there  are  the  errors  inherent 
in  the  physical  process,  due  to  such  factors  as  hydrophone 
variability  or  malfunction,  local  distortions  of  the  sound 
wave,  and  inexactness  of  the  water  column  which  estimates 
the  speed  of  sound  profile.  The  last  two  sources  of  error 
together  make  up  the  variability  that  is  involved  in  the 
time  values  which  are  ultimately  used  in  position  estima- 
tion, and  is  therefore  the  variation  that  is  to  be 
estimated. 

If  it  is  assumed  that  the  variability  indi ced  by  the 
spherical  assumption  is  independent  of  the  data  variability, 
then  the  M.L.S.  variance  estimate  is  the  sum  of  those  two 
variances,    or 

_2  =  —.2  +  2 

Omls  C'sph  U  time 
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Therefore  the  data  variance  can  be  estimated  fcy  first 
estimating  the  variability  induced  by  the  spherical  assump- 
tion, and  then  subtracting  it  from  the  M.L.S.  estimate  of 
the  variation. 

The  M.L.S.  estimate  was  obtained  by  applying  the  M.L.S. 
method  to  the  data  from  a  tracking  run  at  the  Nanocse 
torpedo  tracking  range  on  May  6,  1980.  That  run  involved 
position  estimation  ty  several  different  hydrophone  arrays. 
The  run  made  several  thousand  position  estimates,  733  of 
which  were  at  depths  of  100  feet  or  more  and  involved 
targets  not  more  than  4700  feet  from  the  sensing  array.  The 
depth  limitation  was  imposed  to  avoid  the  excessive  compli- 
cations caused  by  the  radical  changes  in  the  velocity 
profile  above  that  depth  (see  figure  2.  1)  .  The  maximum 
range  limitation  imitates  the  data  validation  procedure  at 
the  tracking  range,  where  positions  farther  than  4700  feet 
from  the  array  are  discarded. 

The  tracking  data  yielded  the  following  range  of  esti- 
mates for  the  standard  deviation  of  the  data  noise,  using 
the  M.L.S.  model. 

MAXIMUM  VALUE  2.89  E-5   sees. 

0.95  QUANTILE  1.18  E-5   sees. 

MEDIAN  VAIDE  3.02  E-6  sees. 

0.05  QUANTILE  4.11  E-7   sees. 

MINIMUM  VALUE  3.13  E-8  sees. 

For  an  overall  estimate  of   the  noise,   the  median  value 

was  used,  so  that: 

rr2  =      (  3.02  E-6  )  2 

u  mis 

The  variability  induced  by  the  spherical  assumption  was 
estimated  by  applying  the  M.L.S.  procedure  to  perfectly 
noiseless  data  in  the  idealized  environment  of  the  single 
array  simulation  of  Chapter  V.  This  was  done  for  targets  at 
ranges    of    1500    to   450C   feet,         at    500    feet    increments,      with 
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1000  randomely  chosen   targets  at  each  range.     The  results 
are  collected  in  Table  XI  . 


TABLE    XI 

Inflation  of   Error   Estimates  Induced 
by  the  Spherical    Model 


Standard   Deviation   Estimates   from   Exact   Tim>     Values 
with   Known   Linear   Velocity    Profile 

RANGE  MINIMUM  Q(.05)  MEDIAN  Q(-95)  MAC  IMUM 

1500  2.72E-10  2.84E-8  2. 08E-7  3.13E-7  3.29E-7 

2000  2.12E-9  5.88E-8  2. 26E-7  3.14E-7  3.30E-7 

2500  3.90E-8  1.09E-7  2.34E-7  3.14E-7  3.3  1E-7 

3000  8.35E-8  1.50E-7  2.37E-7  3.15E-7  3.24E-7 

3500  1.21E-7  1.59E-7  2.37E-7  3.13E-7  3.24E-7 

4000  1.47E-7  1  .65E-7  2. 39E-7  3.16E-7  3.25E-7 

4500  1-46E-7  1 -69E-7  2.42E-7  3.14E-7  3.25E-7 


Table    XI      shows    that      the    median      and    maximun       inflation 

values      are  reasonably     independent   of      target  range.  The 

minimum   values      vary    somewhat/        but   only      for  nnge      values 

below   3000    feet.  This  represents   a    very      stable   situation 

overall.  Therefore   the      inflation      due      to   the      spherical 

assumption    is      estimated   by    the   median      value   at    a      range    of 

3000    feet,    namely: 

rr2    u      =  (    2.37    E-7    )2       . 

V  sph 

Combining  these   two  estimates,   the   varian:  e  estimated 
for  the  timing  data  is 
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2        =         _  2       T  2 

V  time         ^  mis         ^sph 

(3.02  E- 6)2   -   (2.37  2-7)2 
(  3.01  1  E-6  )  2 

As  car,  be  seen,  the  error  induced  by  the  spherical  model 
is  less  than  10%  of  the  H.L.S.  error  estimate.  Therefore 
when  it  is  accounted  for  by  subtraction  from  the  H.L.S. 
variation  estimate,  the  final  variance  estimite  changes 
little. 

The  estimated  value  indicates  a  standard  error  in  the 
6th  decimal  place.  With  a  typical  speed  of  sound  value  of 
4880  feet  per  second,  this  represents  a  position  differen- 
tial of  about 

4880  •  3.011E-6     =     0.015   feet   . 
This  estimate  is  guite  low,   indicating  that  the  time  values 
being  recorded  are  sufficiently  accurate. 

There  is  considerable  opportunity  for  addi:  ional  work 
determining  the  relationship,  if  any,  between  the  time  vari- 
ance and  other  factors  such  as  angles  of  Blevation  and 
azimuth  of  the  target  from  the  array. 

There  is  also  the  problem  that  the  time  / ariation  is 
likely  to  be  array  dependent.  For  example,  consider  figure 
6.1  ,  wherein  the  standard  error  estimates  from  the  actual 
tracking  run  are  plotted  versus  the  range  of  the  target. 
The  plot  does  not  indicate  that  there  is  any  simple  rela- 
tionship between  range  and  error  level.  Howevsr  the  plot 
does  show  a  bunching  pattern.  When  the  error  estimates  are 
plotted  separately  for  each  array,  then  the  bunching  pattern 
becomes  clearly  associated  with  the  individual  arrays. 
Consider  figures  6.2  and  6.3,  where  the  separate  plots  have 
been  made  for  four  different  arrays.  It  is  still  not  clear 
from  these  plots  whether  the  principle  effect  is  due  to  the 
individual  arrays,  or  the  ranges  of  the  targets .  However 
some  level  of  array  dependency  seems  likely,  indicating  a 
need  for  additional  investigation. 
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DATA  FROM  ALL  ARRAYS 
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Figure  6.1    Error  Estimation  Versus  Range  of  Target. 

B.   CHOICE  OF  METHOD 

Clearly  all  indications  are  that  the  plani  r  wavefront 
models,  L.S.  and  M.L.P.  are  not  candidates  for  use  as  posi- 
tion estimators.  Furthermore  the  hybrid  model  L.S.C.  is  an 
interesting  improvement,  but  never  really  performs  well 
enough  compared  to  the  M.L. S.  and  NAVY  models. 

The  original  NAVY  model  is  usually  outperformed  by  the 
adjusted  NAVY_A  method.  However  the  differences  are  not 
always  significant. 

The  spherical  model  M.L.S.,  on  the  other  haid,  consis- 
tently outperformed  all  other  methods  during  the  simulated 
evaluations.  It  would  seem  that  M.L.S.  is  tie  model  of 
choice.  It  does  the  best  job  of  handling  normally  distrib- 
uted errors  in  the   data.    But  that   is  not   tie  strongest 
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DATA  FROM  ARRAY  NO.   14 
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Figure  6.2    Error  Estimation  for  Arrays  7  and  14. 
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DATA  FROM  ARRAY  NO.  56 
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DATA  FROM  ARRAY  NO.  57 
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Figure  6.3    Error  Estimation  for  Arrays  56  and  57 
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argument  for  its  use.  A  more  important,  and  surprising, 
argument  in  its  favor  is  that  when  the  exact,  error  free 
times  are  used,  the  M.L.S.  apparent  position  es timate  will 
usually  produce  the  most  accurate  estimate  of  the  true 
apparent  position.  This  is  the  desired  overall  result,  so 
that  the  sensitive  ray  tracing  procedure  will  be  affected 
minimally  by  apparent  position  estimation  errors.. 

There  are  nevertheless  several  notes  of  caution  which 
should  be  considered  before  embracing  the  M.L.S.  method 
wholeheartedly.  The  first  caution  has  been  stressed  before, 
namely  that  these  conclusions  were  arrived  at  under  the 
idealized  conditions  of  the  simulations  scenarios.  The 
second  caution,  also  previously  stressed,  is  tha:  the  actual 
magnitudes  of  the  differences  between  position  estimates  has 
been  ignored.  It  is  conceivable  that  while  one  method 
always  produces  a  better  estimate  than  another,  the  differ- 
ence between  any  two  position  estimates  is  acceptably  small. 

Lastly  there  is  the  caution  that  the  M.L.S.  model 
involves  a  complicated  iterative  procedure  which  uses 
considerable  computer  time.  It  is  probably  too  slow  a 
procedure  for  use  with  'real  time*  analysis  during  the 
execution  of  tracking  runs.  For  real  time  tracking  the 
NAVY_A  method  currently  in  use  is  probably  preferrable  due 
to  its  simple  computations. 

However,  for  post  run  analysis,  and  also  possibly  for 
calibration  of  the  hydrophone  arrays,  the  M.L.S  method  is 
recommended  as  being  a  more  exact  and  more  robi st  position 
estimator  than  those  methods  currently  in  use. 

C.   RECOMMENDATIONS  FOR  FUTURE  INQUIRY 

Several  recommendations  have  already  been  mi  de  for  work 
needed  to  estimate  the  effect  of  suitable  independent  vari- 
ables on  both  timing  errors  and  the  bias  in  certain  methods. 
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In  addition  there  exist  at  least  two  other  areas  for 
possibly   fruitful  investigation. 

The  first  area  concerns  the  interplay  between  methods. 
Specifically,  the  binomial  comparisons  of  Chapter  V  show 
that  even  the  worst  methods  are  better  than  each  of  the 
other      methods    at      least    part      of      the   time.  Hence    it      is 

possible  that  the  best  method  overall  would  be  a  suitable 
combination  of  methods,  wherein  each  method  is  used  where  it 
is  most  effective.  For  example,  even  though  the  M.L.S. 
method  has  been  indicated  as  the  best  method  for  any 
randomly  selected  position,  it  may  be  consistei tly  outper- 
formed by  another  method  under  certain  circumstances,  such 
as   extremely  high    or   low  elevation   angles. 

The  second  area  for  possible  work  addresses  the  guestion 
of   how      to   next      improve  upon    the      existing   models.  It   is 

herein  suggested  that  the  next  improvement  in  modelling 
would  be  a  method  which  is  based  on  a  linear  velocity 
assumption.  As  figure  2.1  demonstrates,  a  linear  velocity 
profile  is  a  reasonable  approximation  for  most  depths.  This 
would  be  the  next  logical  step  above  the  constant  velocity 
assumption  which  is  associated  with  the  spherical  models. 
Most  of  the  mathematical  basis  for  such  a  model  in  contained 
in  Appendix  A.  Possibly  a  suitable  set  of  equations  could 
be  developed  involving  the  hydrophone  times  and  reflecting 
the  linear  velocity  assumption.  If  so,  the  least  squares  or 
maximum    likelihood    techniques    might    provide   useful    results. 
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APPENDIX    k 
LIHEAE  VELOCITY    PROFILE    THEORY 

All  computations  in  this  study  are  cade  under  the 
assumption  that  the  sound  velocity  is  directly  related  to 
depth 'in  a  linear  manner,  and  is  known  exactly.  Under  these 
circumstances  many  closed  form  results,  not  other  wise  avail- 
able,   can    re  obtained  and   used   in   those   computations. 

Suppose   that   the    velocity    profile    is   given   bf 

V  =  VO  +  \M  *  z 
where  VO  and  V1  are  known  constants,  and  z  is  the  depth 
variable,  measured  down  from  the  water's  surface.  In  this 
case  it  is  known  [Hef.  4]  that  the  path  of  a  soi nd  ray  from 
a  signal  source  to  a  hydrophone  will  have  the  shape  of  an 
arc  of  a  circle.  The  center  of  that  circle  will  be  some- 
where above  the  surface  of  the  water.  The  vertical  place- 
ment of  that  center  is  determined  by  tne  value  of  z  at  which 
the  speed  of  sound  equals  zero  (see  figure  A.1) .  Although 
that  depth  is  negative  and  is  not  really  a  depth  at  all,  it 
nevertheless  has   geometric    meaning. 

Consider  the  vertical  plane  containing  the  circle 
center,  the  sound  source  and  the  hydrophone.  ,et  h  be  the 
variable  which  measures  the  horizontal  position  in  that 
plane.  Let  (h,  z  )=  (a  1,a2)  be  the  position  of  the  hydrophone, 
and  let  (h,  z)  =  (p  1,  p2)  be  the  position  of  the  sound  source. 
Then  C2  given  by  (A- 1)  is  the  Z  coordinate  of  the  circle 
center. 

2  V 

1  (A.1) 

What    must   be   found  is  C1,    the   h    coordinate  of     the  circle 

center,    and   r,    the    radius   of   the   circle.       To   solve    for    these 

values,    the   equation    of    the    circle   is    used,    evaluated    at    the 

the    two   known   locations,    yielding    (A.  2)  . 
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Figure   A.1        Circular  Bay  Paths   for    a    Linear    Velocity   Profile 


2  -  2 


(     ai     "     Cl     >  +         <     *2     -     C2     ) 

(     Pl     "     Cl     )2        +         (     P2     "     C2     )2 


(A. 2) 


The  left  hand  sides  of  the  two  equations  of  (A. 2)  can  be 
equated,  and  solved  for  C1,  leading  to  (A. 3). 

(P0  -  a  ) 


(pl  +  al> 


2(pl  "  al} 


(P   +  a.  -  2C9)     (A.3) 


flhen      this    value      of   Z\       is   substituted      int)     the      first 
equation   of    (A. 2),    then    the    result    is    given    by    (A. 4). 


(a 


i2    ♦ 


-     C     ) 


(A. 4) 


The  circular  arclength  between  the  hydrophone  and  the 
sound  source  is  easily  computed.  Unfortunately  the  velocity 
of  sound  does  not  stay  constant  along  that  arc.  Therefore 
in  order  to  determine  the  amount  of  time  (T)  :  eguired  for 
the  sound  ray  to  travel  the  ray  path,  the  effect  of  the 
velocity  must  be  integrated  along  the  arc.  This  is  dene  by 
(A. 5),  where  S  is  the  arclength  between  the  two  points,  and 
V (s)  is  the  speed  of  sound  as  a  function  of  posL ticn  on  the 
arc. 


d  s 


Vis) 


(A. 5) 


In  (A. 5)  the  sound  source  position  corresponds  to  an 
arclength  of  s  =  0,  and  the  hydrophone  position  correponds 
to  arclength  s  =  S.  In  [  Ref .  4]  this  integral  is  shown  to 
be  equivalent  to  (A. 6). 


da 


cos  u) 


(A. 6) 


In  (A. 6)  A0  is  the  angle  of  elevation  of  the  ray  path  at 
the  sound  source,  and  A1  is  the  elevation  ai gle  at  the 
hydrophone.  The  antideri vative  in  (A. 7)  can  be  used  to 
solve  (A. 6) ,  leading  to  the  ray  path  transit  time  expression 
in    (A. 8)  . 


da 


cos  ( a ) 


In 


1     +    s  in  (  a 
cos ( a ) 


(A. 7) 


In 


"/cos  ( A     ) 
^cos  (a     ) 


1    +    sin(A    )\" 


1    + 


sin  ( Aj ) 


(A. 8) 
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If  the  elevation  angle  at  any  point  along  the  arc  is 
denoted  ty  A,  then  (A. 9)  relates  the  angle  to  thj  derivative 
of  z  with  respect  to  h  along  the  ray  path. 


tan  (a) 


dz 


dh 


(A.S) 


Implicit  differentiation   of  the  equation  of   the  circle 
yields  (A.  10)  as  another  expression  for  the  same  derivative. 

h  -  C, 


dz 

dh 


z  -  C 


Equating  these  two  derivatives  leads  to 
relates  the  elevation  angle  and  the  position 
arc. 

h     -     C, 


tan     (A) 


(A. 10) 

(1.11)        which 
(h,z)       on   the 


(A.  11) 


From    (A.  11)         a    simple      geometric   argument     produces   the 
equations   in    (A.  12). 


z    -    c. 


h     -    C. 


cos      ( A ) 


sin     (A) 


(A.  12) 


First  let  A,  h  and  z  be  equal  to  (Al,a1,a2)  in  (A. 12), 
and  then  let  them  equal  (A0,p1,p2).  Then  substitute  both 
expressions  into  (A. 8).  The  result  is  (A. 13),  the  desired 
expression  for  the  ray  path  transit  time  in  ts rms  of  the 
positions  of  the  sound  source  and  the  hydrophone. 


r  /a 


In 


-  C 


r  +  Pi  -  Cl 
r  +  ai  "  Cl 


(A.  13) 


To  summarize,  if  a  ray  path  is  to  terminate  i  t  the  three 

dimensional  position  (X1#Y1,Z1),   and  the  sound  source  is  at 

(X2/Y2,Z2)/   then   perform  the  following   steps  in   order  to 

calculate  the  exact   ray  path  time  and   elevatioi  angle  that 

would  correspond  to  a  linear  velocity  profile: 
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1.  use  (A.  14)  to  translate  the  three  dimensional  posi- 
tions to  the  two  dimensional  coordinates  of  the 
previously  described  vertical  plane,  with  the  origin 
at  the  water  surface  directly  above  the  end  of  the 
ray  path; 


P1      =   V  (xi  "  x2)2  +  (Yi  "  Y2} 2  (A*  14) 

P2   =   z1       ax   =   0       a2   =   z2 

2.  calculate  C2t    CI    and  r  using  (A.1),  (A.  3)  ,  and  (A.  4); 

3.  calculate  the  exact  ray   path  transit   time  T   using 
(A.  13);  and 

4.  calculate  the  exact   elevation  angle  A  at   the  hydro- 
phone, using  (A. 15). 

A    =    arccos (A.  15) 
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APPENDIX    B 
PAETIA1    DERIVATIVE    FORMULAE    FOR    NEKTON'S    METHOE 

The  central  tool  used  by  Newton's  Method  in  the  develop- 
ment of  the  M.L.S.  model  in  Chapter  IV  is  the  matrix  GP.  It 
is  the  matrix  of  first  order  derivatives  of  the  error 
expressions  g1,  g2,  g3  and  g4.  Those  derivatives  are  set 
forth   in   this  appendix. 

If  X  =  (01/02,03,0),  then  the  error  functions  are 
defined    by     (B.  1)  , 

t.    -    /k~ 
g.(x>         =        c L-.  V    I  i    =    lf2f3         (B.I) 

V*I  M 

T        -     D     M     /     V 

g4(x)         =  u    -    — 

1     -     N 

where      the      values    Ki,      M   and      N    are    the    f unctioi s      given   by 
(B.2). 

2                                                  D2 
K.        =       U        -      (     2UDC./V     )     +     —  (B.2) 

V 


«   -  I 


i-V*? 


3  /T 


M 


3        T  .     -  ,/k~ 
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Recall    that    GP    is  the   matrix    shown    in     (B.3) 


GP    (c   ,C2,C3,U)      = 


u~2 

SC2 
3g4 


dc1       B! 


B^i      Bfi      5^1      B^_ 

3g2 

3g3 

3g4 


SC3 

3C4 

3q2 

a92 

3=3 

3C4 

393 

3?3 

as 

3C4 

3q4 

394 

Be,      Be. 


(B.3) 


Then  given  the  definitions  above,  the  folio/ ing  deriva- 
tives are  the  result  of  straight  forward,  though  tedious, 
differentiation,    and   are   offerred   without    detaila d    proof. 

LEMMA    1 


B  M 

3ci  = 

LEMMA   2 


VK        (T-        K.      )     +     C.T.U     D 

■* ^ i    =    1    2    3       <B'U> 


V     K 


iVKi 


LEMMA    3 


2*u 


C  .T  .       ( 


I10 


3  =  1 


V     K 


3  7^7 


v    U) 


(B.5) 


5n 

3ci 


T  .      U     D 

l 


V     K  .     ,.  /  K 
IV        1 


i     -     1,2,3 


(B.6) 


LEMMA    4 


Bn 
Bu 


I- 


T.      (D    C  .     -     V    U) 


W       VKWK: 


(B-7) 
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FORMULA    1 

The  first  three  diagonal  elements  of  GP  are  the  deriva- 
tives of  each  gi  with  respect  to  Ci  (i=1,2,3)  aid  are  -jiven 
ty     (3.8)  . 

B  M 

(3.8) 


Bg 


3ci 


=  .  i  - 


rT.UDM    -     VK.   (T.      --\/K.      )    -v  „ 


M        K.     ^~ 


FORMULA    2 

The  off  diagonal  elements  in  the  first  three  rows  and 
columns  of  GP  are  the  derivatives  of  each  gi  with  respect  to 
Cj,    and   are   given   by    (B.9). 


a» 


sc3 


i  =  1,2,3 
j  =  1,2,3 
D    ?    i 


(B.9) 


FORMULA    3 

The    derivative    of  each    gi    with      respect   to   U    is    given    by 
(B.  10). 


Hi 


T . M (VU     -     DC . )      +     VK . (T . 
l  111 


.    3  M 

V"^TT  (B.10) 


V     M"     K.    -./    K  . 
i    V        i 


FORMULA    U 

The  first  three  elements  of  the  last  row  of  GP  are  the 
derivatives  of  g4  with  respect  to  each  Ci,  and  i  re  given  by 
(B  .  1  1 )  . 


Bf4  _Bfi 


9  N  cb  M 

(VT„    -    DM)       -       L     -^—        (1-N) 


dc, 


V     (1    -    N) 


(B.  11) 


89 


FCEMU1A  5 

The  last  row,   last  column  of  3P  is  the  derivative  of  gU 
with  respect  to  U,  and  is  given  by  (B.12). 


39. 


h  n 


/3m 


i  - 


(VT  .  -  DM)   -  D  l  ~ 
V   (1  -  N) 2 


t) 


(1  -  N) 


(B.  12) 
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APPENDIX    C 
ONIFORH    SAMPLES    ON    A    TRUNCATED    HEMISPHERE 

The  single  array  simulation  of  Chapter  V  required  a 
random  sample  of  positions  in  space  uniformly  distributed 
over  the  surface  of  a  truncated  hemisphere.  Th2  hemisphere 
is  tc  be  of  radius  r  (3000  feet)  about  the  acoustic  center 
of  a  hydrophonic  array-  The  truncation  of  the  overall 
sphere  is  due  to  the  fact  that  the  upper  portion  of  the 
hemisphere  is  above  the  water  surface,  and  the  lower  half  is 
below    the    sea   bottom. 


dH 


w      = 


r      dH 


=      r   cos(E) 


w     =      r   cos(E)    dH 


Figure    C.1        Hemispherical  Geometry. 

Let  E  be  the  variable  denoting  angles  of  eleration  above 
the  horizontal.  Let  H  be  the  variable  denoting  horizontal 
azimuth  angles  about  the  center  of  the  hemisphere.  Consider 
a   small    piece      of   hemispherical   surface   area      bounded    by    the 
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elevation   angles    E1    and     (E1+dE),         and   by    the   azimuth    angles 
H1    and    (HH-dH)  (see  figure  C.1).         If    H      is   the    horizontal 

width    of    that    piece,    then    H    is    given   by     (C.1). 

w      =      r      dH        =        r   cosCE)    dH  (C.1) 

If  A1  is  the  area  of  that  piece,  then  A1  is  i pproximated 
by  (C-2). 

A1   =    w  dE   =    r  cos(E)  dH  dE  (C.2) 

Now  suppose  that  (  0  <  E1  <  E2  <  Emi x  )  where 
Emax  is  the  elevaticn  angle  of  the  top  of  the  truncated 
hemishpere.  Then  the  ratio  (A1/A2)  of  two  different  areas 
at    elevation  angles    E1    and    E2    is    given  by    (C.  3)  . 

A-j^  r   cos(E    )    dH   dE  cos(E.) 

~~A~J~  r   cos(E2)    dH  dE  cos  (E    )  (c-3) 

If  n1  and  n2  are  to  be  the  (relative)  sampls  sizes  from 
the  two  areas  A1  and  A2  respectively,  then  uniformity  of  the 
sample  requires  that  (C.4)  hold. 

Al  7  nl     =    A2  /  n2  (C.4) 

The   combination    of    (C.U)    and    iC.3)     implies    (:  .  5)  . 

cos (E  ) 
"2        =        nl      cos(E1)  (C5) 

letting  E1  =0,  then  the  relative  sample  size  n2.for 
area  A2  is  given  by  (C.6),  where  n  is  the  relative  sample 
size    at    the   base   of    the   hemisphere. 

n2        =        n   cos(E2)  (C-jfiJ 

Now  the  differential  probability  of  drawii g  a  random 
position   that   has   elevation    angle    S    is   given    by    (C.7),    where 
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N   is    the   total      sample    size    to   be   drawn    from      the    surface    on 
the    sector    bounded    by  the   azimuth    angles    H1    and    (H1+dK). 

fE(E)     "   "T"     =    -Tcos(E)  (c-7) 

Therefore,    if   K    is   defined    to   be    the    constan:     ratio   n/N, 
then    the     corresponding   cumulative   distribution   is      given    by 
(C.  8)  . 


rE 

F    (E)         =         /        K  cos(e)    de  =        K   sin  (E) 


0 


(C.8) 


For  the  highest  elevation  angle,  Emax,  thi  cumulative 
distribution  function  must    equal    one,    so   that    (C.9)    holds. 

F     (E         )         =         1         =        K  sin(E         )  ir    q> 

E     max  max  H-  •  ?) 

As  a  result,  the  constant  K  is  determined  by  (C.10)  . 

K    =    -J 

Sin(Emax)  (C-10> 

Therefore  the  conplete  cumulative  distribute  on  function 
is  given  by  (C.  1  1)  . 

c    r^  sin(E) 

FE(E)         =         "sTali f  (C  11) 

max 

Now  the  inverse  probability  transform  can  be  used.  If  U 
is  a  uniform  random  variable  on  the  interval  (0,1),  then  let 
E  be    given    by     (C.12). 


E     =      arcsin  [ ■ )  ,_     .,„, 

Vsm(E        ))  (C  12) 


max 

Now    choose       an    azimuth    angle      H    randomely      ai d    uniformly 
over    the  interval    (0,27T).       Then    (X,Y,Z)    given  by 

X  =  r  cos  (H)  cos  (E) 
I  =  r  sin(H)  cos  (E) 
Z      =     r   sin(E) 


93 


is  the  corresponding  position  in  spherical  :  oordinat es. 
That  position  is  a  random  position  drawn  from  a  population 
ox  positions  uniformly  distributed  across  the  si rf ace  of  a 
truncated  hemisphere  with  maximum  elevation  angle  Emax. 
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APPENDIX    D 
SINGLE    ARRAY    SIMOIATICN    COMPUTER    PROGR&fl 


C  THIS    FORTRAN    PROGRAM    SELECTS    A    RANDOM    SAMPLE    OP 

C  SIZE    M    OF    3-DIMENSIONAL    POSITIONS    ALL    AT    A    RANGE 

C  OF    F.    FEET    FROM    THE    ACOUSTIC    CENTER    OF    A    HYDROPHONIC 

C  ARRAY.       THE   ARMS    OF    THE    ARRAY    ARE    ASSUMED    TO    BE 

C  ALIGNED    WITH    THE    COORDINATE    SYSTEM    OF    THE    RANGE. 

C  AND    IS   IN    A    POSITION    GIVEN    BY    THE    VECTOR    A. 

C 

C      THE  SPEED  OF  SOUND  PROFILE  IS  ASSUMED  LINE* R, 

C      GIVEN  BY  V  =  VV(1)   +   VV  (2)  *  DEPTH 

C 

C      FOE  EACH  OF  THE  PCSIIONS  SELECTED.  EXACT  Hf DROPHONE 

C      TIMES  ARE  CALCULATED  USING  THE  SUBROUTINE  TCOMP,  AND 

C      RANDOM  ERRORS  ARE  ADDED  TO  THE  EXACT  TIMES.   THE 

C      ERROR  DISTRIBUTION  IS  NORMAL,  WITH  MEAN  ZERO  AND 

C      STANDARD  DEVIATION  SDEV. 

C 

C      THE  TIME  VALUES  ARE  THEN  USED  TO  ESTIMATE  * N  APPARENT 

C      POSITION  BY  TWO  DIFFERENT  METHODS,  USING  APPROPRIATE 

C      SUBROUTINES.   THE  RESULTING  APPARENT  POSITIONS  ARE 

C      THEN  PROCESSED  BY  THE  SUBROUTINE  TRACE  TO  OBTAIN 

C      THE  CORRESPONDING  ACTUAL  POSITION  ESTIMATES.   THE 

C      RESULTING  POSITION  ESTIMATES  ARE  THEN  COMPARED  TO  THE 

C      ORIGINAL  TRUE  PCSITION  TO  SEE  WHICH  METHOD  PRODUCED 

C      THE  MOST  ACCURATE  ESTIMATE. 

C 

C      THIS  IS  DONE  FOR  ALL  POSSIBLE  PAIRS  OF  THE 

C      SIX    METHODS. 

C 

INTEGER  M,I.J,K.N.METH,METH1,METH2,NYES.N¥EST,NN0ISE 
DIMENSION  RELnOOO)  ,RAZ  (1000)  .STUFF  (400  0) 
DOUBLE  PRECISION  P(1000.3)  ,A(3)  ,VV  (2)  ,T  (1000, U) 
DOUBLE  PRECISION  VAR  MOO  0)  ,  TN{1  000  .  H\  ,  NOISE  .  FRACT 
DOUBLE  PRECISION  TEST  ( 1 0 00)  ,PT  ( 1 000 , 3)  . DIF  ( 1 000) 
DOUBLE  PRECISION  DIFT  (1 0 00 5  , AC (3[ , TC  ( 1000) 
DOUBLE  PRECISION  R, Z, SDE V, PI , AZ  (1 000 )  , EL  ( 1 0 00) 
DOUBLE  PRECISION  V,FK ,PHI, PHIMAX, S EED, PEST ( 1 000 , 3) 
C 

C      SET  INITIAL  VALUES  OF  SAMPLE  SIZE,  RANGE, 
C      AND  ERROR  STANDARD  DEVIATION 
C 

M  =  1000 

R  =  3000. DO 

SDEV  =  3.464D-5 


FCRMAT{2X,  '_NOISE  STD.  DEV.   =   «,F15.10) 

6' 9 
SEED  =  931 947. 0D0 


WRITE  (6,  1i  SDEV 
WRITE  (6,99  9 


PI   =    2.D0*DARCOS(0.D0) 
C 

C  SET    ARRAY    POSITION    AND    LINEAR    PROFILE 

C  FCR    SPEED    OF    SCUND. 

C 


A  (1)     =   0.D0 

A  (2)     =   0.D0 

A(31    =    1300.D0 

AC  (1)     =    A  (1)     +     15.0 

AC  (2)    =    A  (2)     ♦     15.0 

AC(3)     =    A(3)     -     15.0 
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c 


c 


c 

c 


VV(1)  =  4840. 7D< 
VV  (2)  =  3314.  D-! 
V    =    VV  (1)  +VV  (2)  : 


=    4840. 7D0 
■5 

.  *A(3) 
C 

C  SET    UPPER     LIMIT   ON    ELEVATION    ANGLE,     IP    NEC2 SSAEY 

C 

IF     (    E    .GT.    A  (3)    )     GO    TO    10 
PHIMAX    =    1.2566D0 
GO    TO    2  0 
10  PHIMAX   =       DAESIN(A  (3) /B) 

C 

C  GENEEATE    RANDOM   VALUES    FOE    ANGLES    OF 

C  ELEVATION     AND    AZIMUTH. 

C 
20  FK    =    1.D0/DSIN  (PHIMAX) 

CALL    GGUBS  (SEED,M,RAZ) 
CALL    GGUBS  (SEED,M, EEL) 
C 

C  CONVERT    ANGLES    TO    3-D    COOEDINATES 

C 

DO    1  1    I    =     1,M 

El  (I)     =    EEL  (I) 

AZJI      =    EAZ    ii 

PHI   =    DAESIN  (EL(I) /FK) 

P(I,1)     =    R*DCOS  (PHI)  *DCOS(AZ  (I)  *PI*2.D0    ^ 

P(I,2J     =    R*LCOS    PHI    *DSIN(AZ  (I)  *PI*2.D0) 


Pjl,3)     =    A  (3)-E*DSIN  (PHI) 
11  CONTINUE 

C 

C  COMPUTE    EXACT     HYDROPHONE    TIMES    UNDEE    THE    LC NEAR 

C  VELOCITY    PEOFILE    ASSUMPTION. 

C 

CALL  TCOMP  (P.  AC,M,VV,TC) 

CALL  TCOMP  (P,A,M,VV,T) 
C 
C      GENEEATE  AND  ADD  EERORS  TO  TIMES. 


NNOISE  =  4*M 
100   CALL  GGNML  (SEED, NNOIS  E, STUFF) 
K  =  1 
DO  111  I  =  1,11 

DO    122    J   =    1,4 

NOISE   =    STOFF(K) 
NOISE    =     NOISE*SDEV 
TN(I,J)     =    T(I,J)+NOISE 
K   =    K    +     1 
122  CONTINUE 

111       CONTINUE 
C 
30  1       FOBMATJ2X,'  (    SAMPLE    SIZE       =       '  ,1)  ,9    )') 


WRITE  (6,36  1) 
-{6,999) 


WRITE 
C 

C  TWO    OUTER    LOOPS    RUN    THROUGH    ALL    PAIES    OF    TIE    SIX 

C  POSSIELS    POSITION    ESTIMATING    METHODS. 


DO    1111    METH1    =1,5 

KMETH    =    METH1    +    1 
DO    1222    METH2    =   KMEIH.6 

IF    (METH1     .  EQ.    METH2)     GO    TO    1222 


200       FOEMATJ2X.  'METHODS    USED    ABE    :     ') 
WRITE  (6,  200) 
METH    =    METH1 
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c 
c 

c 


130 
201 

131 

202 

132 

203 

133 
204 

134 
205 

135 

206 


136 
137 


C 
C 
C 
C 
C 


C 
C 
C 
C 

C 
C 

c 

c 


156 

155 


CALL  SUBROUTINES  TO  IMPLEMENT  THE  METHODS 


IF 


IF 


IF 


IF 


(MET 
CALL 


H  .  NE.  1)  GO  TO  131 
POSNAV  (TN,V,  M, PEST, TEST) 

NAVY, 


FOEMAT (2X. 
WRITE  (6,201) 
GO    TO    150 


uncoere;  TED') 


(METH    . 
CALL    PC 


FOEMAT (2X. 
WRITE  (6,202) 
GO    TO    150 

(METH     .NE.     3 
CA 


NE.     2)    GO    TO    132 
OSNVC     (TN,V,  M, PEST, TEST) 


NAVY,     CORRECTED') 


133 


)     GO   TO 
caLL    POSLS    -(TN,V,M, PEST,  TEST) 
FOEMAT     (18X     'LEAST    SQUARES,     UNCORRECTED') 


WRITE (6,203? 
GO    TO    150 

(METH    .  NE.     4)    GO   TO    134 
CALL    POSLSC     (TN,V,  M#PES'*I 
FORMAT (18X, 'LEAST    SQUARES,     CORRECTED') 
WRITE  (6-204) 
GO    TO     150 


!M.L    POSLSC  .  JTN^V,  M,PESTtTEST) 
IF     (METH    . NE.     5)     GO    TO    135 


CALL    POSMLP     (TN,V, M, PEST, TEST, VAR) 
FOEMAT (18X.  'MAX.    LIKELIHOOD,     PLANAR') 
WRITE  (6,205) 
GO    TO    1 50 

GO  TO  136 


IF  (METH  .NE.  6) 

CALL  POSMLS  (TN ,V, M, PEST, TEST, VAR) 

FORMAT (18x,  'MAX.  LIKELIHOOD,  SPHER 

WRITE (6.206) 

GO    TO    150 
WRITE  (6-137)     METH1.METH2 
FORMAT(2X, 'CANT    FIND 


ICAL' 


GO  TO  TOO 6 


METHODS  ',14, •  AND/OR  ',14) 


THE  APPARENT  POSITION  ESTIMATES  ARE  RELATIVE  TO 
THE  ARRAY  CENTER,  AND  MUST  BE  TRANSLATED  TO 
1HE  TRACKING  RANGE  COORDINATE  SYSTEM. 

'150   DO  144  I  =  1,M 


144       CONTINUE 


PEST(I,1)  =  PEST  (1,1)  ♦  AM 
PEST'1,2)  =  PEST(I/2[  +  A(2 
PESTiI,3)     =    A(3)     -    PEST  (1,3 


CONVERT    APPARENT    POSITIONS    TO    ACTUAL    ESTIMATES    BY 
RAY    TRACING    PROCEDURE. 

CALL    TRACE     (PEST,PT,  TEST, A , M, VV) 

CALCULATE    DIFFERENCE    BETWEEN    THE    1ST    POSITION 
ESTIMATE    AND    THE   TRUE    POSITION. 

IF     (    METH     .NE.     MET 31     )     GO    TO    160 
DO    155    I    =    1,M 
DIF  (I)    =    0.D0 

DIFTfl)     =    DABS  (TEST  (I) -TC  (1,4)) 
DO    156    J    =    1,3 

DIF(I)     =    DIF(I)  +  (PT  (I,J)-P  (I,  J)  )  **> 
CONTINUE 
CONTINUE 
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METH    =    METH2 
GO    TO     1  30 
C 

C  CALCULATE    DIFFERENCE    BETWEEN    THE    2ND    POSITION 

C  ESTIMATE     AND    TEE    TRUE    POSITION,     AND    COMPAR2     I' 

C  TO    THE    1ST    POSITION    DIFFERENCE. 

C 
160       NYES    =   0 
NYEST    =    0 
DO    166      I    =    1,M 
DO    167    J  =    1,3 

DIF(I)     =    DIF(I)  -  (PT  (I,J}-P(I,J))  **2 
167  CONTINUE 


DIFT(I)     =    DIFT(I)     -    DABS  (TEST(I)-TC  (1,4)  ) 
'Da FT (I)    .  GT.     0.nn 
YEST   =    NYEST*  +    V 


165  IF     (    DIF(I)     .GT.    0 . DO    )     GO    TO    166 

NYES    =    NYES    +    1 

166  CONTINUE 

FRACT    =     (DFLOAT  (NYES) /DFLOAT(M)) 
300       FORMAT (2X,  f    FRACTION    OF    TIME    METHOD    #1    IS    :  LOSER') 
303       FOEMAT(2X,'  IN       POSITION  :       ',F7.3) 

302       FORMATJ2X,'  IN      TIME  .       \F7.3) 

WRITE  (6,300) 

WRITE  (6,30  3)     FRACT 
=     (DFLOAT 

WRITE  (6,30  2)     FRACT 

WRITE  (6,999 

WRITE  (6,998 

WRITE  (6,  999 
C 

1222    CONTINUE 
1111    CONTINUE 
C 

998  F0EMAT(1X,'  =    =    =    =    =    =    =    =    =    =    =    =    =    ==    =    =    =    ^) 

999  FORMAT(2X,'  ') 
1000       SIOP 

END 


??ACT    =    ]DFLOAT_(NYEST) /DFLOAT  (M)) 
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APPENDIX  E 
DOUBLE  ARRAY  SIMULATION  COHPaTER  PROGR&M 


C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c 

c 


c 
c 
c 
c 


THIS    FORTRAN    PROGRAM    SELECTS    A    RANDOM    SAMP- E    OF 
SIZE    M    OF    3-DIMENSIONAL    POSITIONS    IN    A    BOX    RUNNING 
CROSSWAYS    BETWEEN    TWO    HYDROPHONIC    ARRAYS.        THE    BOX 
HAS    DIMENSIONS    GIVEN    BY    THE    VECTOR    BOX.       THE    ARRAYS 
ARE    SEPERATED    EY    7500     FEET.    AND    HAVE    COORDINATES 
GIVEN    BY    THE    VECTORS    A1    AND    A2.       THE    ARMS    3 F    THE 
ARRAYS    ARE    ALIGNED    WITH    THE    RANGE    COORDINATE    SYSTEM. 

THE    SPEED    OF    SOUND    PROFILE    IS    ASSUMED    TO    BE    LINEAR, 
AND    IS   GIVEN    BY  V      =       VV(1)     +    VV(2)*DEPTH 

FOR    EACH    OF   THE   POSITIONS    SELECTED,    EXACT    TIME    VALUES 
FOR    SCUND     WAVE    ARRIVAL    AT    THE    HYDROPHONES    OF    EACH 
ARRAY    ARE    COMPUTED    BY    THE    SUBROUTINE    TCOMP.        THEN 
RANDOM   ERRORS    ARE    ADDED    TO    ALL    TIMES.       THE    ERROR 
DISTRIBUTION    IS   NORMAL,     WITH    MEAN    ZERO    AND    STANDARD 
DEVIATION    SDEV. 

THE    TIME    VALUES    ARE    THEN    USED    TO    ESTIMATE    IPPAEENT 
POSITIONS    BY    TWO    DIFFERENT    METHODS.     USING    THE 
APPROPRIATE    SUBROUTINES.        EACH    METHOD    PRODUCES    TWO 
DIFFERENT    ESTIMATES,     ONE    FOR    EACH    ARRAY.       THE    APPARENT 
POSITIONS     ARE    THEN    TRANSLATED    TO    ACTUAL    POSITION 
ESTIMATES    BY    THE    SUBROUTINE    TRACE.       THEN    THE    LENGTH 
OF    THE   DIFFERENCE    VECTOR    FOR    EACH    METHOD    IS    COMPUTED. 
THE    TWO    DIFFERENCE    LENGTHS    ARE    COMPARED    TO    SEE    WHICH 
METHOD   PRODUCES    POSITION    PAIRS    IN    CLOSEST    AGREEMENT. 


THIS    COMPARISON    IS 
THE    SIX    METHODS. 


DONE    FOE    ALL    POSSIBLE    PURS    OF 


INTEGER    M/I/ J,K,METH,  METH1 , ME TH2, N YES, NNOIS  E,N ARRAY 
DIMENSION    RP  (1000)  , STUFF  (4  0  00)  ,  NAM  1  (6)  ,  NAM 2  (6) 


DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 
DOUBLE 


PRECISION 
PRECISION 
PRECISION 
PRECISION 
PRECISION 


A1  (3)  ,A2(3)  ,T1  (10  00,4)  ,T2(  1000,4) 
VARf  100  0)  .NOISE. FR ACT, XYZ. DEPTH 
TEST  (10  0OJ  ,?  1(1000  ,3)  ,P2(f  0  00,4) 
SDEV,  SEED,  BOX  (3)  ,  PEST  ( 1000  ,  3)  ,  A  (3) 
T  (10  00,4)  ,V,VV(2)  ,DIF(1000) 


DATA  NAM1 

DATA  NAM  1 

DATA  NAM 2 

DATA  NAM2 


)  ,NAM1  | 

[2) 

,  NAM1 | 

[3) 

)  ,NAM1 

5) 

,  NAMU 

6 

,NAM2i 

2 

,  NAM2 

[3 

) ,NAM2 

5) 

,  NAM 2 

6 

/'NAVY' , 'NAVY', 'L.S.» 
/•L.S.  '  ,  «M.  L.  •  ,  •  M.L.  • 
/•  •f'CORR','  ' 

/'CORR',  'PLNR',   'SPHR' 


INITIALIZE    VALUES    FOR    SAMPLE    SIZE, 
AND    ERROR     STANDARD    DEVIATION. 


RANGE,  DEPTH 


M  =  1000 
SDEV  =  3.464D-5 
DEPTH  =  1300. DO 
FOEMAT(2X,'  NOISE 
WRITE  (6,  1[  SDEV 
WRITE  6.999) 
SEED  =  934  7.  6D0 


STD.  DEV. 


,F15.  10) 
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c 

c 


30 


c 
c 
c 


c 
c 

c 


22 
11 


c 
c 
c 
c 


33 


45 

44 


56 
55 


SET    UP    ARRAY    POSITIONS,    AND    SOUND    VELOCITY    PROFILE 

A1  (1)     =   -3750. DO 
A1    2      =    0.D0 
A1  (3)     =    DEPTH 
A2  (1)    =    3750. DO 
A2  (2)     =    0.D0 
A2  (3)    =    DEPTH 

=    4840. 7D0 
•5 

(1)  +  VV  (2)  *DEPTH 


az  io  i     -    jjx.  ri  n 
VV  (1)     =    4840. 7D( 
VV    2)     =    33  14.  D-i 
V   =    VV/1)  +VV  12)  •■ 


FORMAT  (2Xf  • 
WRITE  (6,30)  M 
WRITE  (6,  999) 


(  SAMPLE  SIZE   =   ',!> 


)  f) 


DRAW  UNDERWATER  BOX  FOR  RANDOM  POSITIONS 
BOX(1)  =  1000. DO 


BOX(1)  =  1000. DO 
BOXi'2)  =  5000. DO 
BOX  (3)  =   DEPTH 


GENERATE  RANDOM  POSITIONS  IN  BOX 

DO  1  1  J  =  1,2 

CALL  GGUBS  (SEED,M,RP) 

DO  22  I  =  1  ,  M 
XYZ  =  RP  (I) 
P1(I,J)  =  (XYZ-  0.5D0)  *BOX  (J) 

CONTINUE 
CONTINUE 

CALL  GGUBS  (SEED, M,RP) 
DO  33  I  =  T,M 

XYZ  =  RP(I) 

P1  (I,  3)  =  XYZ*BOX(3) 
CONTINUE 

COMPUTE  EXACT  HYDROPHONE  ARRIVAL  TIMES,  AND  ADD 
RANDOM  ERRORS  TO  TIME  VALUES. 

NNOISE  =  M*4 

CALL  TCOMP  (P1 ,A1,M,VV ,T1) 

CALL  GGNML  (SEED, NNOI SE, STUFF) 

K    =    1 

DO    44    I    =     1,M 

DO    45    J    =    1.4 

NOISE    =     STUFF(K) 
NOISE    =    NOISE*SDEV 

+ 


T1  (I.J)     =   T1  (I,  J) 
K    =    K+1 


NOISE 


CONTINUE 
CONTINUE 


CALL    TCOMP  (P1  ,A2,M,VV,T2) 
CALL    GGNML     (S EED, NNOI SE, STUFF) 
K    =    1 
DO    55    I    =    1,M 

DO    56    J    =     1,4 

NOISE   =    STUFF(K) 
NOISE   =     NCISE*SDEV 
T2  (I, J)     =   T2  (I,  J)     +    NOISE 
K   =     K+1 
CONTINUE 
CONTINUE 
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c 
c 


10 

20 


C 
C 

c 


125 


C 
C 


67 
66 


130 
131 
132 
133 
134 
135 


136 
137 


C 
C 

c 
c 
c 


144 


C 
C 

c 


SET    OF    LOOPS    TC    BUN    THHOUGH    ALL    PAIRS    OF    M2THODS. 

DO    1111    KETH1    =1,5 

KMETH    =    METH1    +    1 
DO    1222    METH2    =    KMETH. 6 

IF     (METH1     .  EQ.    METH2)     GO    TO    1222 


FORMAT (2X, 'METHODS  COMPARED  ARE  1 
FOBMAT(2X     '  AND       2 

WRITE  (6,  16)  NAM1  (METH  1)  ,NAM2(METH1 
WRITE  (6,  20j  NAM1  (METH2J  ,NAM2(METH2 
WRITE  (6,999) 

MEIH    =    METH1 
NAERAY   =     1 

SET    UP    ARRAY    EEING    USED 


' ,A4,2X,A4) 
',A4,2X,A4) 


(1}     =   A1  (1) 
2      =   A1     2 
3)     =   A1     3) 


A 
A 
A 
DO  66    I    =  "1,M 

DO    67    J    =     1,4 

T(I.  J)     =    T1  (I,  J) 

CONTINUE 
CONTINUE 

CALL    SUBROUTINES    TO   PERFORM    ESTIMATION    METHODS. 

IF 


(METH    .NE.     1)    GO    TO    131 
CALL    POSNAV     (T , V, M ,PEST, TESI ) 


IF 


IF 


II 


IF 


GO    TO     150 
IF     (METH    .  NE.     2)    GO   1 0    132 

CALL    POSNVC     (T,V,M  , PEST, TEST) 

GO    TO    150 

(METH    -NE.     3)    GO    TO    133 

CALL    POSLS      (T,V,M,  PEST,  TEST) 

GO    TO     150 

(METH    .  NE.     4)    GO    TO    134 

CALL    POSLSC     (T,V,M, PEST, TEST) 

GO    TO     150 

(METH    *NE.     5)    GO    TO    135 

CALL    POSMLP     (T,V,M,  PEST,  TESI,  VAR) 

GO    TO     1  50 

H    .  NE.     6)    GO   TO    136 
POSMLS     (T,V,M  ,  PEST,  TEST,  VAR) 

GO    TO    150 
WRITS  (6.  137)      METH1,METH2 

FORMATfSX.'CANT    FIND    METHODS    ',14,'     AND/OR     »,I4) 
GO    TO    1000 


(MET! 
CALL 


150       IF     (    NARRAY    . EQ.    2   )     GO    TO    152 

APPARENT    POSITION    ESTIMATES    ARE    IN    LOCAL    AJRAY 
COORDINATES,     AND    MUST    BE    TRANSLATED    TO    TRACKING 
RANGE    SYSTEM    COORDINATES. 


DO    144   I    =    1,M 

PEST(I,  1)     =    PEST  (1,1 
PESTfl,  2)     =    PEST(I 
PESTJI,3)     =    A1  (3) 
CONTINUE 


,1)     +    A1  (1) 
,2(    +    A1     2 
-    PEST  (1,3 


CORRECT    APPARENT   POSITIONS    BY    RAY    TRACING. 
CALL    TRACE    (P  EST,  P  1,  T  EST  ,A  ,  M,  VV) 
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c 
c 


78 
77 


C 

c 
c 
c 


152 


145 


C 

c 
c 


c 
c 

c 


c 
c 
c 
c 


156 
155 


GO    ON    TO    SECOND    ARRAY 
NARRAY   =    2 


AM)     =   A  2(1) 

A  (2)     =   A2    2) 
AJ3)     =   J2    3 


'h 


160 

167 

166 


300 
310 


C 
1222 
11  1  1 

C 
99  8 

999 
1000 


I    =     1,M 

DO    78    J    =     1,4 

Tfl-  J)     =    12  (I,  J) 

CONTINUE 
CONTINUE 
GO    TO    130 

TRANSLATE    2ND    ARRAY    APPARENT    POSITIONS    TO    1ANGE 
COORDINATE    SYSTEM,    AND   THEN    RAY    TRACE. 

DO    145   I    =    1,M 

PEST(I,1)     =    PEST(I,1)     +    A2  (1 

PEST  i'I,  2)     =    PEST  (1,2)     +    A2  (2 

PESTi'1,3)     =    A2(3)     -    PEST  (1,3 
CONTINUE 

CAIL    TRACE     (PEST, P2, T EST , A, M, VV ) 

COMPUTE    DIFFERENCE    VECTORS    FOR     1ST    METHOD 

IF     (    MEIH    .NE.     METH1     )     GO    TO    160 
DO    155    I    =    1,M 
DIF(I)    =    0.D0 
DO    15  6    J    =    1,3 

DIF(I)     =    DIF(I)+ (P1  (I,J)-P2(I,J))  **2 
CONTINUE 
CONTINUE 

GO    CN    TO    SECOND   METHOD 

METH    =    METH2 
NARRAY    =    1 
GO    TO    125 

COMPUTE    DIFFERENCE    VECTORS    FOR    2ND    METHOD    *ND 
COMPARE    TO    DIFFERENCES    FOR     1ST    METHOD. 

NYES    =    0 
DO    166      I    =    1 ,M 
DO    167    J    =     1,3 

DIF(I)     =    DIF(I)-(P1(I,J)-P2(I,J))**2 
ITINUE 


CONTI! 

brkf 


IF    {    DIF(I)     .GT.    O.DO    )     GO    TO    166 
NYES    =    NYES    +    1 
CONTINUE 


METHOD    #1    ') 

:     '  ,F7.3   = 


WRITE  (6,99  9) 

FRACT    =     fDFLOAT(NYES)  /DFLOAT(M 

FOBMAT^X,'    FRACTION    OF    TIME 

FORMAT]2X.«  IS    BETTER    IS 

WRITE  (6,300) 

WRITE  (6,310)     FRACT 

WRITE  (6,999) 

WRITE    6,  998) 

WRITE  (6,  99  9) 

CONTINUE 
CONTINUE 

FORMATMX,^    =    =    =    ===    =    =    =    ===    =    =    ===•) 

FORMAT  (2  X,'  f) 

STOP 
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APPENDIX    F 
COMPUTER    SUBROUTINES    FOR    TIME    CALCULATION    AMD    RAITRACING 


C    SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS5  ssssssss 

c 

SUEROUTINE    TCOWF    (P,A,M,VV,T) 
C 

C  THIS    FORTRAN    SUEROUTINE    COMPUTES    THE    EXACT    TIME 

C  REQUIRED    FOR    A    SOUND    WAVE    TO    TRAVEL    FROM    ITS 

C  SOURCE    (VECTOR    F)     TO    THE    FOUR    HYDROPHONES    ON    AN 

C  ARRAY    WHOSE    ACOUSTIC    CENTER    IS    SPECIFIED    3Y 

C  VECTOR   A. 
C 

C  THE    METHODOLOGY    USED    IS    SET    FORTH    IN    APPENDIX    A 

C  OF    THE   IHSSIS.       THE    BASIC    ASSUMPTION    IS    ONE    OF 

C  A    LINEAR    VELOCITY    PROFILE,     WHOSE    COEFFICIENTS    ARE 

C  GIVEN    BY    THE    VECTOR    VV. 


C 


INTEGER    M.I.J 


DOUBLE  PRECISION   P  { 1  000 ,3)  ,  A  (3)  ,  VV  1 2)  .  T  ( 13  0  0,  4) 
DOUBLE  PRECISION   AA ( 4,3) , P 1 , P2 , R, CI , A2, C2 


C 

C      SET  UP  COORDINATES  OF  FOUR  HYDROPHONES  ON  THE  ARRAY 
C 

DO    11    I    =     1,4 

DO    22    J    =    1  ,3 

AA  (I.J)     =  -15.  DO    +    A  (J) 
22  CONTINUE 

AA(I,3)     =     15. DO    +    A  (3) 
11  CONTINUE 

AA  (1,  1)     =    15. DO   +    A(1) 
AA  (2,2)     =    15.  DO    +    A(21 
AA(3;3)     =    -15. DO    +    A  (3) 
C 

C      CALCULATE  QUANTITIES  C1,  C2,  R,  AND  THE  Til  ES  T(*,4) 
C 

C2  =  -1.D0*(VV  (1)/VV(2)  ) 
C 

C  IOCP    THROUGH    THE    M    SOURCE    LOCATIONS 

DO    33    I    =     1,M 
P2    =   P  (1,3) 
C  LOOP    THROUGH   THE    FOUR    HYDROPHONES    FOR    Ei CH    SOURCE 

DO    44    J    =     1,4 

P1    =    DSQRT(  (P(I  ,1)  -AA  (J,  1)  )  **2 

♦    (Pa,2)-AA(J.2))  **2) 
C1    =    P1**2+P2**2-AA{J,3)  **2 

+  2.D0*C2*(AA(J,  3)-P2) 
C1    =    C1/(2.D0*P  1) 

~T  (C1**2  +  (AA  (J,3)-C2)  **2) 
((AA(J,J)-C2J/(P2-C2)  ) 

*  ((R-C1+P1)/(R-CI  )) 
(DLOG(T(I,J)))/VV(2) 


R    =    DSQR 
T(I,J)     = 

44 
33 
C 

T(I,J)     = 
CONTINUE 

CONTINUE 

RETURN 
END 
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C    SS3SSSSSSSSSSSSSSSSSSSSSSSSS35SSSSSSSSS3SSSSSSSSSSSSSSSS 

c 

SUBROUTINE       TRACE     (P,  PT  ,  T  ,  A,  M,  V  V) 

c 

C  THIS  FORTRAN  SUBROUTINE  TAKES  AN  APPARENT  POSITION 

C  ESTIMATE  P  AND  CONVERTS  IT  TO  AN  ACTUAL  POSITION 

C  ESTIMATE  PT  BY  RAY  TRACING.   THE  RAY  TRACING  ASSUMES 

C  A  LINEAR  VELOCITY  PROFILE,  WITH  COEFFICIENTS  IN  THE 

C  VECTOR  VV.   OTHER  INPUTS  ARE  THE  ARRAY  LOCATION 

C  VECTOR  A  AND  THE  RAY  TRANSIT  TIME  TO  THE  ACOUSTIC 

C  CENTER.  T.   THE  RAY  TRACING  LAYERS  ARE  DEL  FEET  THICK 

C  WITH  THE  FIRST  (BOTTOM)  LAYER  BEING  THE  FIRST  DEL 

C  FEET  IMMEDIATELY  ABOVE  THE  ARRAY. 

C 

C 


INTEGER   M,  I 

DOUBLE  PRECISION  T  (10  00)  ,P(  10  00,3)  ,PT  (1000,  3)  ,A  (3) 
DOUBLE  PRECISION  DEL, H,S A, CA, DEPTH , V , V 1 , VV ( 2) 
DOUBLE  PRECISION  DS,DT,DZ,DH 


DEL    =    25. DO 
C 

C  LOOP    THROUGH     THE    M    APPARENT    POSITIONS 

C 

DO    11    I    =     1,M 
C  INITIALIZE     INCREMENTAL    VALUES    FOR    EACH    POSITION 

H    =   DSQRT  I  (P  (I,  1)  -AM))  **2+  (P  (1,2)  -A  (2)  i   **2) 

CA    =    H/DSQRT(H**2  +  (P  (I,3)-A(3)  )  **2) 

DT    =    0. DO 

DZ    =    0.D0 

DH    =    0. DO 

DEPTH    =    A (3)    -    DEL/2. DO 

V  =   VV  (1)  +VV  (2)  *DEPTH 
C 

C  INNER    LOOP    :     PROCESS    DATA    UPWARDS,    LAYER    Bt     LAYER, 

C  UNTIL    RAY    TRANSIT    TIME    IS    EXHAUSTED 

C 

10  SA    =    DSQRT  (  1.D0-CA**2) 
DZ    =    DZ    +    DEI 

DS    =    DEL/SA 
DT    =    DT    +    DS/V 
DH    =    DH    ♦    DS*CA 
C 

IF     (    DT    .GT.    T  (I)     )     GO    TO    20 
DEPTH    =    DEPTH   -    DEL 
V1    =    VV{1)     +    VV(2)  *DEPTH 
CA    =    CA*(V1/V) 

V  =    V1 
GO    TO     10 

C 

C  ADJUST   FOR    OVERSHOOTING    IN    LAST    LAYER 

C 
20  DS    =    V*  (DT-T  (I)  ) 

DH    =    DH    -    CA*DS 
C 

C  ADJUST   APPARENT    POSITION    TO    GET     ACTUAL    POSC  TION 

C 

PT(I,1)     =     (P(I,1)-A(1j)*(DH/H)+A(1) 

c 

11  CONTINUE 
FETURN 
END 
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APPENDIX    G 
COMPUTER   SUBROUTINES    FOR    POSITION    ESTIMATION    METHODS 


C    SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS5SSSSS 

c 

SUBROUTINE       POSNAV     (T,V,M,P,NT) 

c 

INTEGER       M,    I 
C 

C  THIS    FORTRAN    SUBROUTINE    IMPLEMENTS    THE    ORIGINAL 

C  UNADJUSTED    NAVY   APPARENT    POSITION    ESTIMATION    METHOD. 

C 

C  INPUT    ARE    THE    HYDROPHONE    TIMES    T    FOR    M    SOU:*  D    SOURCE 

C  POSITIONS,     AND    VELOCITY    V    AT    THE    ARRAY. 

C 

C  OUTPUT   ARE    M    APPARENT    POSITION    ESTIMATES    P    ALONG    WITH 

C  THE    CORRESPONDING    M    RAY    TRANSIT    TIMES    NT.        ALL    THE 

C  POSITIONS     ARE    REFERENCED    TO    THE    ACOUSTIC    C2NTER,     WITH 

C  THE    Z    COMPONENT    BEING    MEASURED    UPWARDS    FROM    THE    ARRAY. 

C 

DOUBLE   PRECISION    T  (10  00/  4)  ,  P  ( 10  00  ,  3)  ,  NT  (1 0)  0) 
DOUBLE   PRECISION    V,TC ,D, DISC, NU MER, DENOM 
C 

D    =   30. DO 
DO    11    I    =     1,M 
NUMER    =    0.DO 
DENOM    =    0.D0 
TC    =    T(I,4) 
DO    22    J   =    1,3 

P(I,J)     =    V*V*(TC*TC-T  (I, J)  **2)  /(D*2.D0) 
NUMER    =    NUMER    +    P(I,J)**2 
DENOM   =    DENOM    +     (P(I,J)     +    15.D0)**2 
22  CONTINUE 

NT  (I)     =    TC*DSQRT(NUMER/DENO«) 
T1  CONTINUE 

FETURN 
END 
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c 


THIS    FORTRAN     SUBROUTINE    IMPLEMENTS    THE    ADJUSTED    NAVY 
METHOD    NAVY_A    FOR    ESTIMATING    APPARENT    POSITIONS. 


C    SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS5  sssssssssss 

c 

SUEROUTINE    POSNVC     (T,V,M,P,NI) 

c 
c 
c 

c 

C  INPUT    ARE    THE    HYDROPHONES    TIMES    I    FOR    M    SOJ ND    SOURCE 

C  LOCATIONS,    AND    THE    VELOCITY    V    AT    THE    ARRAY. 

C 

C  OUTPUT   ARE    THE    M    APPARENT    POSITION    ESTIMATES    P,    AND 

C  THE    CORRESPONDING    M    RAY    TRANSIT    TIMES    NT.        ALL 

C  POSITIONS    ARE    REFERENCED    TO    THE    ACOUSTIC    CENTER    OF 

C  THE    ARRAY,    WITH   THE    Z    COMPONENT    MEASURE    UPWARDS 

C  FROM    THE    ARRAY. 


INIiiGEh    M,I 

DOUBLE   PRECISION      T  (1  000,4)  ,P(1  000-3)  -NT  (13  001 

DOUBLE   PRECISION      V,  D  ,  NUMER  ,DENOM,  TC  ,  DCC 


D    =    30. DO 
DO    11    I    =     1,M 
DCC   =    O.DO 
TC    =    T(I,4) 
DO    22    J   =    1,3 

P(I,J)     =    15. DO    +     (V*V*(TC*TC-T  (I,  J)*"  2)  /(D*2.D0)  ) 
DCC    =    DCC    +    P(I,J)**2 
22  CONTINUE 

NUMER  =  O.DO 
DENOM  =  O.DO 
DO    Hf|    J    =     1,3 

PJIcJ)     =     (P  (I,J)  *V*IC)/DSQRT  (DCC) 
DENOM    =    DENOM    +    P(I,J)**2 
P(I,J)     =    P(I,J)     -    15. DO 
NUMER    =     NUMER    ♦    P(I,J)**2 
44  CONTINUE 


11  CONTI 

RETURN 
END 


NT  (I)     =    TC*DSQRT(NUMEfi/DENOM) 
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c  ssssssssssssssssssssssssssssssssssssssssssssssssssssssssss 
c 

SUBROUTINE    POSLS     (T, V , M, P, 1ST) 

c 

C  THIS    FORTRAN     SUBROUTINE    IMPLEMENTS    THE    LEAST    SQUARES 

C  PLANAR    METHOD    I.S. 

C 

C  INPUT  ARE  THE  HYDROPHONE  TIMES  T  FOR  M  SOUND  SOURCE 

C  POSITIONS,  AND  THE  VELOCITY  V  AT  THE  ARRAY. 

C 

C  INPUTS   AND    OUTPUTS    ARE    THE    SAME    AS    FOR    THE 

C  SUBROUTINE    POSNAV. 


C 


c 


INTEGER  M,I,J 

DOUBLE  PRECISION  T (10 00, 4) ,P ( 10 00 , 3) , LST ( 1) 00) 

DOUBLE  PRECISION  V,DISC,TC 

DO  11  I  =  1.H 
DISC  =  0.D0 
DO    22    J    =    1,3 

P(I,J)     =    T(I,4)     -T{I,J) 
DISC    =    DISC    +    P (I, J) **2 
22  CONTINUE 

IST(I)     =     (T  (I,1)+I(I,2)+T(I,3)-T{I,4))/2.  DO 
DO    33   J    =    1 ,  3 

P(IrJ)    =      (V*LST  (I)  *P  (I,  J)  )  /DSQRT(DISC) 
33  CONTINUE 

CONTINUE 

SETUBN 
END 
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C  SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS53SSSSSSSS5  sssssssssss 

c 

SUBROUTINE    POSISC     (1,  V#H#P#1ST) 

c 

C  THIS    FORTRAN    SUBROUTINE    IMPLEMENTS    THE    BIAS    ADJUSTED 

C  LEAST    SQUARES     PLANAR     METHOD    L.S.C.     FOR    ESTIMATING 

C  APPARENT    POSITIONS. 

C 

C  INPUTS   AND    OUTPUTS    ARE    THE    SAME    AS    FOR    THE 

C  SUBROUTINE    POSLS. 

C 

INTEGER    I,M,J,K.L 

DOUBLE   PRECISION    1(10  00,4)  ,P(  1000,3)  ,151(1)  0  0) 

DOUBLE   PRECISION    V , A ( 4, 3) , D  (4 )  , DISC, TD, R IS QR , R2SQR, E (3) 
C 
C  SET    UP   HYDROPHONE    POSITIONS 

DO    44    J    =1 ,3 

DO    55    I    =    1,4 

A  (I,  J)     =    -15. DO 
CONTINUE 


A(J.J)      =    15. DO 

-our 


55 

44         CONTINUE 
C 

C  CALCULATE    ORIGINAL   L.  S.     SOLUTION 

DO  11  I  =  1,  M 
DISC  =  0.D0 
DO    22   J    =    1,3 

FjI.J)    =    1(1,4)     -TCI«J1 
DISC    =   DISC    +    E>  (I,J)**2 
22  CONTINUE 

LSTtt)     =    (T(I,1)+T(I,2)*T(I,3)-T(I,4))/2.  DO 
DO    33    J    =    1,3 

P(I,J)     =     (V*LST(I)  *P  (I, J))  /DSQRT  (DISC) 
33  CONTINUE 

R1SQR_=    P(I,  1)**2+P  (1,2)  **2  +  P  (1,3)  **2 

DO    89    K    =    1,3 

D(4)     =    D(4)     +     (P  (I,K)     -    A(4,K))**2 
89  CONTINUE 

D]4)    =    DSQRT  (D(4)  ) 

DISC   =    0.D0 

TD    =   0.D0 

DO    77    J    =    1  , 3 

D(J)     =    0.D0 


DO    88    K    =     1,3 

D(J) 
CONTINUE 
D(J)     =    DSQRT  (D(J)  ) 
DISC    =    DISC    +     (D  (4)  -D  (J)  )  **2 


D(J)     =    D(J)     ♦     (P(I,K)     -    A(J,K))**2 

TINU" 


TD    =    TD    +    D(J) 

77  CONTINUE 

C 

C  CALCULATE    BIAS    VECTOR    E 

DISC   =     DSQRT  (DISC) 
DO    66    J    =    1,3 

E(J)     =     (  (1D-D  (4)  )  *  (D  (4)-D(J)  )/(DISC*2.  ,D0)  )  -P  (I, J) 
66  CONTINUE 

R2SQR    =    0.D0 
DO    99    J    =1,3 

PJI, J)     =    P(I,J)     -    E  (J) 
R2SQR    =    R2SQR    +    P(I,J)**2 
99  CONTINUE 

C 
C  ADJUST    ORIGINAL    RAY    TRANSIT    TIME 

1ST  (I)     =    LST  (I)  *(DS  QRT  (R2SQR/R1SQR)  ) 
11  CONTINUE 

RETURN 
END 
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C    SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS3SSSSSSSSSSSSSSSSSSSSS 

c 

SUBROUTINE       POSMLP     (T  ,  V,  M,  P,  MT,  VAR) 

c 

C  THIS    FORTRAN     SUBROUTINE    IMPLEMENTS    THE    MAXC  MOM 

C  IIKELIHOOD    PLANAR    METHOD    M.L.P.     FOR    ESTIMATING 

C  APPARENT    POSITIONS. 

C 

C      INPUIS  ARE  THE  HYDROPHONE  TIMES  T  FOR  M  SOURCE 

C      POSITIONS,  AND  THE  VELOCITY  OF  SOUND  V  AT  THE  ARRAY. 

C 

C      OUTPUTS  ARE  THE  APPARENT  POSITION  ESTIMATES  P  AND 

C      RAY  TRANSIT  TIMES  MT  FOR  EACH  OF  THE  M  POSITIONS. 

C      ALSO  OUTPUT  IS  A  VECTOR  VAR  OF  M  ESTIMATES  OF  THE 

C      TIMING  ERROR  VARIANCE.   POSITIONS  ARE  ALL  REFERENCED 

C      TO  THE  ACOUSTIC  CENTER  OF  THE  ARRAY,  WITH  THE  Z 

C      COMPONENT  MEASORED  UPWARDS  FROM  THE  ARRAY. 

C 

INTEGER       M,    I,     J 

DOUBLE   PRECISION    T  (1  0  00,  4)  ,  P  (  10  00 ,  3)  ,  MT  (1  0)  0) 
DOUBLE   PRECISION    V, TC ,D, DIFF, NUMER,DENOM,T3 L, C  (3) 
DOUBLE   PRECISION   V  AR  (  1 00  0)  ,  CO  (3)  ,  X  (3) 
C 

D    =   30. DO 
TOL    =    1.D-4 
C 

C  INITIALIZE    THE    DIRECTION    COSINE    ESTIMATE    Bf 

C  USING    THE    LEAST    SQUARES    SOLUTION. 

C 

C  LOOP    THROUGH    THE    M    POSITIONS 

C 

DO    11    I    =     1,M 
TC    =    T(I,4) 

DENOM    =     (TC-T(I,1)  )**2+  (TC-T  (1,2)  )  **2 

♦(TC-T    1,3       **2 
DO    12    J    =    1,3 

C(J)     =      (TC    -    T (I, J)) /DSQRT (DENOM) 
12  CONTINUE 

C 

ITER    =    0 
22  DENOM    =    0.D0 

ITER    =    ITER    '+    1 
TC    =    0. DO 
C 

C  USE   ITERATION    FORMULAE    TO    DEVELOPE    COSINES 

C  AND   TIME    VALUES. 


C 


DO    33    J    =    1,3 


DENOM    =    DENOM    +    T(I,J)*C(J) 
TC    =    TC    +   T  (I,  J) 
33  CONTINUE 

TC   =     ((TC+T  (1,4) )  +  D*(C[1)+C(2)  +C(3))/V)  /4.  DO 

DENOM    =    DENCM    -    TC  *  (C  (1)  +  C  (2  )  +C  (3)  ) 
DO    66    J    =    1,3 

C(J)    =     JT  (I,J)-TC)/DENOM 

XJJ)     =    DABS  (CO  (  J) -C(J)) 

66  CONTINUE 

C 

C  CHECK    ITERATION    TOLERANCES 

C 

DIFF   =    DMAX1  (X(1)  ,X(2)  ,X  (3)  ) 
IF     I    ITER     .LT.     100)     GO    TO    70 

60  FORMAT     (2X, ' ITE RAT  IONS    EXCEED    ',16) 

61  FORMAT     (2X,,IN    POSITION    NO.  • ,  1 6) 
WRITE(6,60)     ITER 

WRITE(6,61)     I 
GO    TO    11 
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70  IF     (    TOL    .11.    DIFF    )     GO    TO    22 

C 

NUMER    =    0.D0 
DEHOH    =    0.D0 
C 

C  CONVERT    COSINES    TO    POSITIONS,     AND    TRANS,  ATE    TO 

C  ACOUSTIC   CENTER    REFERENCED    COORDINATES. 

C 

DO    44    J    =     1,3 

Pfl-  J)     =    V*C  (J)  *TC 
DENOfl    =    DENOM    +    P(I,J)**2 
P(I,J)     =    P(I,J)     -15.D0 
NUMER    =     NUMER    +    P(I,J)**2 
44  CONTINUE 

MT  (I)     =    TC*DSQRT(NUMER/DENOM) 
C 

C  ESTIMATE    VARIANCE 

C 

VAR  (I)  =  0.  DO 
DO  55  J  =  1,3 

VAR  (I)  =  VAR  (I)   +  (T  (I,  J)  -TC  +  D*C(J)  /I)   **2 
55       CONTINUE 

VAR(I)  =  VAR(I)/4.  DO 
11    CONTINUE 
RETURN 
END 
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c  ssssssssssssssssssssssssssssssssssssssssssssssssssssssssss 
c 

SUEROUTINE      POSMLS     (T  IMS  ,  V  ,M,  P,  'ALT,  VARML) 

c 

C  THIS    FORTRAN     SUEROUTINE    IMPLEMENTS    THE    MAXC  MUM 

C  LIKELIHOOD    SPHERICAL    METHOD    M.L.S.    FOR    ESTIMATING 

C  APPARENT    POSITIONS. 

C 

C  INPUTS   AND    OUTPUTS    ARE    THE    SAME    AS    FOR    THE 

C  SUEROUTINE    POSMIP. 

C 

INTEGER       M.    I-    J,    ISR,    II.     ITER1,     ITER2 
DOUBLE   PRECISION    TIME  { 1  0  00  ,  4)  ,  P  (1  0  00  ,  3)  ,  VAX  ML  { 100  0) 
DOUBLE   PRECISION    ML!  (  1000)  ,  SDEV  M  000)  ,  K  (3)  ,     DSDC(3) 
DOUBLE   PRECISION    V,R , D ,TOL 1 . IOL 2. EPS . L, S, DIF F 
DOUBLE   PRECISION    T  (3)  ,  TC  ,TAU,  TK  (3)  ,  VK  (3  )  ,  C  {  3)  ,  C1  (  3) 
DOUBLE   PRECISION    XJ3     , DENOM, NUM ER, DL DC  7 3)  . D LDT r DSDT 
DOUBLE   PRECISION    GP  (4  ,  4)  ,GPI  (4,  4)  ,  WKARE  A  (2  8  )  ,  GC  (4) 
DOUBLE   PRECISION    GT,  X  1  (4)  ,  X2  (4)  ,  A  (4)  ,  B  (  4)  ,  T  1 
DOUBLE   PRECISION    F1,F2 

C 

C         SET    GLOEAL    VALUES 


C 


R    =    (3.D0-DSQRT(5.D0)  )/2.D0 


D    =   30. DO 

TOL1    =    1.D-5 

TOI2    =    1.D-5 
C 
C         START    OUTERMOST    LOOP     (ONE    FOR    EACH    SOURCE    POSITION) 


C 


DO    5    II    =     1,M 
DENOM    =    0. DO 
TAU   =    TIMEjII,4) 
EPS    =    1.D-3 
IC    =   TAU 
DO     111    J    =    1,3 


T(J)     =    TIME(II,J) 

CJJ)     =     (  (V**2)*(TC**2-T[J)  **2)/(D*2.D)  ))  +D/2.D0 

DENOM    =    DENOM    +    C(J)**2 


C(J)  **2 

111  CONTINUE 


CJJ)     =    C  (5)/DSQRT  (DENOM) 
NUE 


C 

DO    122    J    =    1#3 

122  CONTI 

C 

C         START    MIDDLE    LOOP       (CREATE    DERIVATIVE    MATRIX    G?(,)     ) 
C 

ITER1    =    0 
10  11    =   TAU 

ITEE1    =    ITER  1    ♦    1 

IF     (ITEE1    .  LT.    50    )     GO    TO    300 

310  FORMAT(2X, 'EXCESSIVE    NEWTON    ITERS.) 

311  FORMAT(2X,*  IN    POSITION       *,I6) 
WRITE'(b#3  10) 

WHITE  (6,311)  II 

GO  TO  5 
300     L  =  0.D0 
S  =  0.D0 
DLDT  =  0.D0 
DSDT  =  0.D0 
C 

DO    11    I    =    1.3 

C1  (I)     =    C  (I) 

K(I)     =    TAU**2+  (D/V)  **2-(2.  D0*D*TAU*C  (C  )  ) /V 

VK(I)     =    V*K(I)*DSQRT(K(I)  ) 

IK]l)     =    T  (I)-DSQRT(K(in 

DLDC(I)     =     (  (V*K  (I)  *TK  (I)  )  + 


DLDC(I)  =  (  (V*K  (I)  *TK  (1)  ) +Ci 
DSDC(I)  =  IT  (I)  *D*TAU)  /VK(I) 
DLDT    =    DLDT+  (C(I)*T  (I)  *  (D*C 


C(I)  *I  (I)  *D*IAU)/VK(I) 
I) 
C'(I )  *T  (I)  *  (D*"C  (I)  -V*TAU)  )  /VK  (I) 
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c 


c 


c 


DSDT    =    DSDT+T  (I)*  (D*C  (I)  -V*TAU) /VK  (I) 
L=    L+(C  (I)*TK(I)  )/DSQRI  (K  ' 
£   =    S+TK  (ij/DSQRT  (K  (I)  ) 


\  "    L+iC  (I)*TK_(I).)^DSCJfiI  (K(I)  ) 

CONTINUE 


DO    22    I    =    1,3 

DO    33    J    =     1,3 


GP  (I,J)  =  (TK(I)  *DLDC  (J)  )  /  (L*L *DSQRT(  K  (I)  )  ) 
33  CONTINUE 

GP  (4  , 1)  =DSDC  (I)  *  (V*TC-D*L) 

GP(4,I    =  (GP(4,I)  -D*DLDC  (I *(1  .DO-S)  ) 

/(V*(1.  DO-S)  **2) 
GP(I,I)  =  (1*T(I)  *D*TAU)  -DLDC(I)  *K  (I)  *V*TK  (I) 

") *L*L)t 

*K(I)  *C  K  (I)  *DLDT 


GP(I,IJ     =     1.D0    -     (GP(I,I)/(VK  (I)  *L*L)i 
GP(I,4)  =  (T  (I)  *L*  (V*TAU-D*C  (I)  )  ) 

+V*K(I)  *( 


GP(I,4)     =    GP(I,4)    /     (VK(I)*L*L) 
22  CONTINUE 

GP  (4,4)  =  1.D0-((DSDT*(V*TC-D*L))  -  D* DID T* M . DO-S)  ) 

/(V*(t.  DO-S)  **2) 
C 

C         INVERT    DERIVATIVE   MATRIX 
C 

CALL   3AUSS3  (4,0, GP, GPI,I2R,  4) 
C 

C  CALCULATE    INITIAL    NEWTON    SEARCH    SOLUTION 

C 
40  DO    44    I    =    1 ,3 

GC(I)     =    C  (I)-TK(I)/  (DSQRT(K(I)  )  *L) 
44  CONTINUE 

GC(4)     =    TAU- (V*TC-D*L)/(V*  (1.D0-S)  ) 
B(4)    =    TAU 
DO    55    I    =    1.3 

B  (4)     =    B  [4>-GPI(4  ,1)  *GC  (I) 

DO    56~J    =     1,4 

B(I)    =     E(I)-GPI(I,J)*GC(J) 
56  CONTINUE 

55  CONTINUE 

B(4)    =    B  (4)     -   GPI(4  ,4)  *GC  (4) 
C 
C  PREPARE   FOR    GOLDEN    SECTION    SEARCH 


EPS    =    DMAX1  (M.D-6)  ,  (EPS/1.  D1)) 
DO    66    I    =    1,3 


CONTINUE 


&j)«MgrB*(B(1)-1(1» 


A(4)    =    TAU 

XT  (4)     =    A  (4)  +E*  (B(4)  -A  (4)  ) 
TAU   =    X1  (4) 
C 

C  START    INNERMOST    LOOP     (    COMMENCE    GOLDEN    SECTC ON    SEARCH 

C  TO    IMPROVE    THE    INITIAL    NEWTON    SOLUTION    ) 


ITER2    =    0 

70  CALL   OBJFCN     (F1 , L, K ,TK,S , X1 ,D , V,T, TC) 

ITER2    =    ITER2    +    1 
IF     (ITER  2    .LI.    50    )     GO    TO    4  00 
410  FORMAT^X, 'EXCESSIVE    G.S.S.     ITERS.,    P)  S    #',I6) 

WRITE  (6,4  10)     II 
GO    TO    5 
400  DIFF   =    0.D0 

DO    77    I    =    1,3 


DIFF    =    DIFF    +     (A  (I)-3(I)  )**2 
X2  (I)     =    A  (I)     ♦    B  (I)     -    X' 
C(I)     =    X2    I) 

NTINUE 


B  (I)     -    X1  (I) 
77  CO 
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X2  (4)     =    A(4)     +    B  (4)     -    XI  (4) 
TAU   =    X2  (4) 

DIFF   =    DSO.RT(DIFF    +  (A  (  4)  -B  (  4)  )  **2) 


C 

C     CHECK  TOLERANCES  -  STOP  G.S.S.  IF  TIGHT  ENOJ GH 


C 


IF     (   EPS    .GT.    DIFF    )    GO    TO    100 

CALL   OBJFCN (F2,L,K,  TK, S, 12, D, V, T, TC) 


C 

C  CHOOSE    IMPROVING   DIRECTION 


C 


IF     (   F1     .LT.     F2    )     30    TO    90 
DO    88    I    =    1,3 


88  CO 


Am    =  x*ui) 

XT(I)     =    X2(I) 
C(I)     =    X1  (I) 

NTINUr 


A(4)     =    X1j4) 
XI  (4)     =     X2(4) 

GO    TO    7  0 
C 
90  DO    99    I    =    1,3 


BJI)     =    X^  (I) 

X1(I)     =    A  (Ij  +B(I)-X1  (I) 

C(I)     =    X1     I) 

99  CONTINUE 
B(4)    =    X2(4) 

XI  (4)     =    A(4)  +B{4)  -X  1  (4) 

GO    TO    70 
C 

C  CHECK    TOLERANCES    FOR     NEWTON    SEARCH    ITERATI) NS, 

C  AND    PREP    FOR    NEXT    NEWTON    ITERATION    IF    NECESSARY 

C 

100  DO     144    I    =     1,3 

X(I)     =    DAES  (C  (I)  -C1  (I)  ) 
144  CONTINUE 

DIFF   =    DMAX1  U(1)  ,X  (2)  ,X  (3)  ) 
IF    {   TOL1    .LT.    DIFF    )     GO    TO    10 
DIFF   =    DABS  (TAU-T1) 
IF     (   TOL2    .LI.    DIFF    )     GO    TO    10 
C 

C  DONE    WITH       II-TH   SET    OF    TIMES    AND    POSITIONS.         MAKE 

C  ESTIMATES,     AND    GO    ON    TO    NEXT    POSITION    TO    BE    ESTIMATED 

C 

VARML(II)    =    0.D0 
DENOM    =    0. DO 
NUMER    =    0. DO 
DO    155    J    =    1,3 

VARML(II)     =    VARML(II)     +    TK  (J)  *TK  (J) 
P(II.J)     =    V*TAU*C(J) 
DENOM    =    DENOM    +    P(II,J)**2 
P(II,J)     =    V*TAU*C  (J)-D/2.D0 
NUMER    =    NUMER    +    P(II,J)**2 
155  CONTINUE 

VARMLUI)     =      (VARML(II)  +  (TC-TAU)  **2)/4.D0 
SDEVJII)     =     DSQRT  (VARML  (II)) 
MLT  (II)     =    TAC*DSQRT  (NUMER/DENOM) 
5  CONTINUE 

C  WRITE(7.900)      (SDEV(I),I    =    1,M) 

C900  FORMAT  (5E15.  5) 

EETURN 
END 
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C  SSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSSS5SSSSS5SSS5  sssssssssss 

c 

SUEROUTINE  OBJICN  (F,  I,  K,TK  ,S  ,  X  ,  D  ,  V,  T,  TC) 

c 

C  THIS  FORTRAN  SUBROUTINE  IS  CALLED  BY  THE  SJBROUTINE 

C  POSMLS.   IT  CALCULATES  NEW  VALUES  FOR  SEVERAL 

C  VARIAELES  AS  WELL  AS  A  FUNCTIONAL  VALUE  WHICH  IS  THE 

C  DECISION  FACTOR  DETERMINING  THE  APPROPRIATE 

C  IMPROVING  DIRECTION  FOR  THE  GOLDEN  SECTION  SEARCH. 


C 


INTEGER  I 

DOUBLE  PRECISION  F,L,  K  (3)  ,TK  (3)  ,  S,  X  ( 4)  ,  D,  V,  TC  ,T  (3) 

S  =  0.D0 
L  =  0.D0 
DO  11  I  =  1,3 

K  (I)  =  XY4)  **2+  (D/V)  **2-((2.  D0*D*X  (4)  *X(  I)  ) /V) 

TK  (I)  =  T(I)-DSQRI  (K(I)  ) 

S  =  S  +  TK  (I)/DSQRT(K  (I)  ) 

L  =  L  +  (X  (I)*TX(I))/DSQ£I  (K  (I)) 
11    CONTINUE 
F  =  0.D0 
DO  22  I  =  1,  3 

F  =    F+  (X  (I)  -TK  (I)  /  (DSQRT  (K  (I)  )  *L)  )  **2 
22    CONTINUE 

F  =  F+(X  (4)-  (V*TC-D*L)/(V*  (1.D0-S)  )  )  **2 

RETURN 

END 
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